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TO ALL WHOM IT MAY CONCERN: 

Be it known that we, Rafael Yuste, Vikram S. Kumar, Robert C. Froemke, and 
Paul Czkwianianc, citizens of Spain, the United States, the United States, and Poland, 
residing in the States of New York, Massachusetts, California and New York, whose post 
office addresses are: 601 West 1 1 3th Street Apt. 5G New York, NY 1 0025; 390 
Commonwealth Avenue, Apt. 605, Boston, MA 02215; 1977 Pleasant Valley #6, 
Oakland, CA 9461 1; and 126 East 56th St., 19 Floor, New York, NY 10022, respectively, 
have invented an improvement in: 

METHOD AND SYSTEM FOR ANALYZING FINANCIAL MARKET DATA 

of which the following is a 

SPECFICATION 

RELATED APPLICATION 
[0001] This application claims priority from U.S. provisional application No. 

60/245,132 filed on November 2, 2000, which is incorporated by reference herein in its 

entirety. 

BACKGROUND OF INVENTION 

[0002] The present invention relates to analyzing and interpreting datasets of financial 
market information. Examples of such datasets include closing price information for 
multiple financial instruments over time. As used herein, financial instrument means any 
commodity, security, instrument or contract traded on an open or closed market or 
exchange including stocks, bonds, options, future contracts, promissory notes and 
currencies. 
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[0003] It is often desirable to understand the relationship of various events occurring 
within a financial market information dataset. For example, share prices for various 
stocks may rise or fall with certain cohesiveness. It is desirable to determine which, if 
any, group of stocks ever exhibited correlated behavior (i.e. share prices rise or fall at the 
same time at least once in the period of observation), regularly exhibited correlated 
behavior (i.e. share prices rise or fall together on multiple occasions over the period of 
observation), and which stock, if any, consistently rises or falls before or after another 
stock rises or falls. It would also be advantageous to know the statistical significance of 
the relationships between the various events. In other words, whether the correlation 
among the various events is stronger than would be expected from random activity. 

SUMMARY OF THE INVENTION 
[0004] These and other advantages are achieved by the present invention which in one 
respect provides a method for analyzing a financial market dataset and for detecting 
relationships between various events reflected in the dataset. 

[0005] In an exemplary embodiment, a method is presented for analyzing a financial 
market data array with a first dimension and a second dimension. The array is examined 
to detect events of interest, and those events of interest are stored in an event array having 
the same dimensions as the financial market data array, but the data in each element of 
the event array is binary. The financial market data array or the event array is then 
analyzed to determine relationships between the events of interest and correspondingly, 
relationships between the financial instruments corresponding to the financial market 
data. 
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[0006] In an additional exemplary embodiment, analyzing includes plotting a portion 
or all of the data in the first simplified array to allow visual examination of the 
relationships between the activities of interest. In another exemplary embodiment, the 
analysis step involves detecting events of interest that are coactive and determining 
whether the number of coactive events is statistically significant. This embodiment may 
include detecting all such coactive events (i.e. instances where events where events occur 
in at least two financial instruments simultaneously), detecting instances where many 
financial instruments are coactive simultaneously, or detecting instances where two or 
more financial instruments are each active in a certain temporal relationship with respect 
to one another (also referred to as coactivity). 

[0007] In a further exemplary embodiment, the data analysis involves calculating a 
correlation coefficient between two financial instruments based on how often the 
financial instruments are coactive relative to how often the first financial instrument is 
active. Representations of all such financial instruments are displayed with lines between 
representations of the financial instrument having a thickness proportional to the 
correlation coefficient between the two financial instruments. 

[0008] Another exemplary embodiment includes plotting a cross-correlogram or 
histogram of events of interest in a particular financial instrument with respect to events 
of interest in another financial instrument, so that the histogram will reveal the number of 
times an event of interest in the first financial instrument occurs a certain number of 
locations away from an event of interest in the second financial instrument. The cross- 
correlogram can be plotted with respect to only one financial instrument, thus showing 
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how many times an event of interest occurs before or after the occurrence of another 
event of interest in the same financial instrument. 

[0009] Yet another exemplary embodiment includes displaying a time series "movie" 
showing activity occurring in one or more financial instrument relative to activity in a 
selected financial instrument. This "movie" is referred to herein as a spike triggered 
average. In this embodiment, a number of frames before and after events occurring in the 
selected financial instrument is chosen. A movie having the number of frames chosen is 
then displayed, with icons displayed for each non-selected financial instrument that was 
active within the chosen number of frames before or after activity occurring in the 
selected financial instrument. A parameter of the icon for each non-selected financial 
instrument, such as the color of the icon, is varied in each frame of the movie to 
correspond to the frequency that non-selected financial instrument is active and the 
corresponding number of frames before or after events occurring in the selected financial 
instrument. 

[0010] Other exemplary embodiments include performing Hidden Markov Modeling 
on the event array to determine a hidden Markov state sequence and displaying a cross- 
correlogram between events of interest occurring in one region of interest while that 
region is in one of the detected Markov states and performing a singular value 
decomposition on the financial market data array. 

[0011] In another aspect of the present invention there is provided a system for 
carrying out the foregoing method. 
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BRIEF DISCRETION OF THE DRAWINGS 
[0012] For a more complete understanding of the present invention, reference is made 

to the following detailed description of exemplary embodiments with reference to the 

accompanying drawings in which: 

Fig. 1 illustrates a flow diagram of a method in accordance with the 
present invention; 

Fig. 2 illustrates a visual plot generated in accordance with the method of 

Fig. 1; 

Fig. 3 illustrates an example of a data structure useful in the method of 

0 Fig. 1; 

p Fig. 4 illustrates a flow diagram of a method of analyzing data useful in 

l"U the method of Fig. 1; 

|t| Fig. 5 illustrates a visual plot generated in accordance with the method of 

* Kg- 1; 

M» Fig. 6 illustrates a cross-correlogram generated in accordance with the 

p method of Fig. 1; 

Jjif Fig. 7 illustrates a correlation map generated in accordance with the 

M method of Fig. 1; 

Fig. 8 illustrates an exemplary format for displaying analysis results useful 
with the method of Fig. 1; 

Fig. 9 illustrates another exemplary format for displaying analysis results 
useful with the method of Fig. 1; 

Fig. 10 illustrates yet another exemplary format for displaying analysis 
results useful in the present invention; and 

Fig. 11 illustrates yet another exemplary format for displaying analysis 
results useful in the present invention. 
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DESCRIPTION OF THE PREFERRED EMBODIMENT 
[0013] Referring to Fig. 1, there is shown a flow diagram representing an exemplary 

method for analyzing data pertaining to financial instruments in accordance with the 

present invention. For purposes of this description, the financial instrument data is 

arranged in an input array corresponding to a time series of daily closing prices for 

various publicly traded stocks. Thus, the data array is a two dimensional array, with one 

dimension (indexed by a first dimensional index) corresponding to the different stocks 

and the other dimension (indexed by a second dimensional index) corresponding to the 

dates the closing prices were observed. The format of this input data array will be 

discussed further herein with reference to Fig. 3. It will be understood that the present 

invention is not limited to the particular data described. For example, the input data 

could correspond to any parameter of any type of financial instrument sampled at any 

frequency. For example, rather than including closing price data, the input data array 

could consist of price/earning ratios, market capitalization or trading volume of the 

various stocks over time. Alternatively, the data could consist of closing quoted prices 

for a commodity, such a electricity, available for delivery at a certain geographic 

location. Moreover, rather than consisting of daily closing prices, the data could consist 

of prices observed at the expiration of any other temporal period, such as every five 

minutes, or every month. Numerous other potential input data sets will be apparent to 

one of ordinary skill in the art. 

[0014] In the exemplary embodiment, performance of the method is assisted by a 
general purpose computer with a processor adapted to operate the MAC-OS operating 
system and to interpret program code written in Interactive Data Language ("IDL") 
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version 5.1 or later, developed by Research Systems, Inc. The IDL program code of the 
exemplary embodiment is appended hereto as Appendices A, B and C described further 
herein. Other operating systems and programming languages could be used to perform 
the steps of the exemplary embodiment without departing from the scope of the 
invention, and the modifications necessary to make such a change will be apparent to one 
of ordinary skill in the art. 

[0015] In step 101, events of interest in the input financial data array are detected. To 
further understand this step in the exemplary embodiment, reference is made to Fig. 3 
where an example of an input data array 300 is shown. Data array 300 is a two 
dimensional array input data having multiple rows 322, 324 . . . 326 and multiple columns 
321, 323 ... 325. Each one of the rows 322, 324 ... 326 corresponds to a particular 
financial instrument, such as a particular stock. Thus, all data within a single row consists 
of observations corresponding to the same stock. Although only three rows are shown in 
Fig. 3, it will be understood that any number of rows could be present, the number of 
rows corresponding to the number of stocks under analysis. Each one of the columns 
321, 323 ... 325 corresponds to a particular time period, such as a particular day on 
which the observation was made. Thus, all data within a single column consists of 
observations occurring during the same day. Although only three columns are shown in 
Fig. 3, it will be understood that any number of columns could be present, the number of 
columns corresponding to the number of observations made. Each data element, 301, 
303, 305, 307, 309, 311, 313, 315, 317 corresponds to a particular observation. For 
example, data element 309 corresponds to the observation of the stock corresponding to 



FILE NO. AP33682 (70050.1603) 

PATENT 

EXPRESS MAIL LABEL NO: EF321680474US 
row 324 made during the period corresponding to column 323. Thus, data element 309 
may contain the closing price of stock A observed on day X . In that scenario, data 
element 307 (which is in the same row as element 309) would contain the closing price of 
stock A observed during the period corresponding to column 321 and data element 315 
(which is in the same column as element 309) would contain the closing price of the 
stock corresponding to row 326 observed on day X. 

[0016] To assist in comparing the observations of different financial instruments 
trading at different prices, the data in input matrix 300 may be modified to contain 
percent change observations rather than actual closing price observations. For example, 
the closing price information for the stock associated with each row 322, 324 . . . 326 of 
input data could be modified to contain percent change rather than absolute closing prices 
as follows. Beginning with the data element in the second column 323, the difference in 
closing price from the observation in first column 321 to the observation in second 
column 323 is calculated. The resulting difference is then divided by the closing price 
observation in the first column 321. The resulting value is stored in the data element in 
the second column 323. The process is repeated until the final column 325 is reached. 
Each element in the first column of data (i.e. data elements 301, 307 ... 313) is then set to 
zero. In this fashion, each data element will represent the percent change in closing price 
from the previous observation, rather than containing raw closing price data. 

[0017] Returning now to Fig. 1, in step 101 the events of interest in the input data array 
300 are detected. In one exemplary embodiment an event of interest is detected by 
calculating a statistical mean and standard deviation for all data elements corresponding 
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to a particular stock. Thus, where the input data is contained in the array 300, a mean and 
standard deviation is calculated for all data in each row of the simplified array. An event 
is then detected where the data element value exceeds the mean for all data in the row by 
a predetermined number of standard deviations. If activity were defined by a drop in 
value rather than an increase in value, the event could be detected by examining the data 
values in a financial instrument for an entry where the data element value is less than the 
mean for all data in the row by a predetermined number of standard deviations. The 
number of standard deviations may be entered by a user before the calculations are 
preformed, or a default number may be used, such as two or three. In this fashion, the 
method will detect those instances in time where the closing price is much higher than the 
average closing price, thus suggesting an event of interest has occurred. 

[0018] In another exemplary embodiment, an event is detected by looking for a data 
value that exceeds a previous data values corresponding to the same stock instrument by 
a threshold amount. Thus, for example, if the closing price stored in data element 309 
exceeded the closing price stored in data element 307 by a certain percentage, an event is 
said to have occurred at the time corresponding to data element 307. Again, if an event 
were indicated by a drop in value rather than an increase, the detection step would 
involve looking for a stock price that is less than previous stock price of the same stock 
by the threshold amount. The threshold amount can be specified by a user before the 
calculations are performed, or a default number can be used, such as five percent. The 
detection can occur over many time periods, for example, the closing price of a particular 
stock on day six could be compared to the stock's closing price on day one to see if an 
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increase beyond the threshold amount has occurred over that period. This would be 
useful to detect events that occur gradually over time rather than relatively 
instantaneously. 

[0019] In step 103, the results of detection step 101 are stored in an event array. For 
this purpose, the event array is identical to the input array illustrated in Fig. 3; however, 
the data stored in the event array is binary rather than closing price values or percent 
changes. Thus, the entries in the event array would be 1 or 0 (or yes or no), 
corresponding to whether an event of interest occurred in the corresponding stock at the 
corresponding time. 

[0020] In step 105, the stored data is analyzed. In one exemplary embodiment, the data 
is analyzed to determine whether various stocks are correlated (i.e. whether they are 
coactive), the strength of those correlations (i.e. how often they are coactive relative to 
how many times each stock or one of the stocks is active), how significant the 
correlations are (i.e. whether the correlation is stronger than would be expected if from a 
random data set) and the behavior of the entire observed stock population. 

[0021] In the exemplary embodiment, the data is analyzed by plotting at least a portion 
of the data contained in the input data array 300. For example, stock price for one stock 
can be plotted over time. Stock prices for all observed stocks could also be plotted over 
time, either in separate plot windows or superimposed on the same plot window in either 
two or three dimensions. Additionally, the closing prices for all stocks could be averaged 
and plotted over time to show global behavior of the observed stocks. Fig. 2 illustrates 
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one possible plot of stock closing price over time, expressed as percent change as 
previously described. 

[0022] In another exemplary embodiment illustrated in Fig. 5, the data is analyzed by 
plotting at least a portion of the data contained in the event array. As shown, a plot of 
events over time may be presented for one or multiple stocks in the input data set. For 
example, events occurring in three stocks are shown plotted versus time in Fig. 5. Events 
for each stock are plotted on separate horizontal axes 501, 503 . . . 505. The vertical lines 
507, 509, 511 represent events occurring at respective times in the corresponding stock. 

[0023] In yet another exemplary embodiment illustrated in Fig. 4, the data in the 
financial data array is analyzed to determine the number of coactive events in the dataset 
and the statistical significance of those events. In step 401, a random distribution of 
stock price activity is generated. The random data is generated by shifting the data in 
each row of the input data array by a random amount. In step 403, the number of 
coactive events in the random dataset is counted. This process is repeated numerous times 
to generate a random distribution. The number of random trials may be set by the user or 
a default number of random trials may be conducted, such as 1000. 

[0024] Counting coactive events for this purpose means counting all instances where 
two stocks are coactive. Coactive events for this purpose means events of interest that 
occurred in two stocks at the same time, or within a specified number of time intervals 
from each other. Thus, if the specified number of time intervals is one, then if a event 
occurred in the stock corresponding to row 322 at the time corresponding to column 321 
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(i.e. data element 301) and an event occurred in the stock corresponding to row 324 at the 
time corresponding to column 323 (i.e. data element 309), those events would be 
considered coactive. The time interval may be specified by a user before coactive events 
are counted, or may be a default setting such as two time intervals. 

[0025] Once the random trials have been completed and a random distribution of 
coactive events generated, the actual number of coactive events in the data is calculated 
in step 405 using the same counting methodology was used to count coactive events in 
the random trials. The actual number of coactive events is then superimposed on a plot 
of the random distribution. The statistical significance of the coactive events is 
determined in step 407 by calculating the area under the distribution curve to the right of 
the number of actual coactive events in the data. This result, termed the "p-value" 
represents the probability that the number of detected coactive events in the actual data is 



[0026] In a further exemplary embodiment, a random distribution of activity is 



generated as previously described, except the only coactive events that are counted in 
steps 403 and 405 are those where a predetermined number of stocks are coactive. The 
predetermined amount of coactive stocks may be specified by a user or a predetermined 
default value such as four may be used. Additionally, it may be specified whether 
exactly that many coactive events must be present or at least that many coactive events 
must be present to be considered a coactive event for counting. Thus, the embodiment 
allows instances of multiple simultaneously active stocks (rather than simply two 
simultaneously active stocks) to be counted and the statistical significance of that number 




produced by a random activity. 
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to be reported. In this exemplary embodiment, the random distribution and actual 
number of coactive events are plotted. The statistical significance of the actual number of 
coactive events is calculated using the formula: C ra nd/N ran d where C ra nd is the number of 
random trials that resulted in more coactive matches than the actual data set and N ran d is 
the total number of random trials used to generate the random distribution, and is 
reported to a user. Additionally, a chart may be drawn showing all observed stocks with 
line segments connecting those stocks that were coactive, such as the chart described 
herein with reference to Fig. 7. 

[0027] In a still further exemplary embodiment, a random distribution of stock activity 
is generated as previously described except the only coactive events that are counted in 
steps 403 and 405 are those where at least two stocks are active a predetermined number 
times throughout the dataset. The number of times the two or more stocks must be active 
can be specified by a user or a default number such as two may be used. In this 
exemplary embodiment, the random distribution and actual number of coactive events are 
plotted. The statistical significance of the actual number of coactive events is calculated 
using the formula: C ra nd/N ra nd where C ran d is the number of random trials that resulted in 
more coactive matches than the actual data set and N rand is the total number of random 
trials used to generate the random distribution, and is reported to a user. Additionally, a 
chart may be displayed showing all observed stocks with line segments connecting those 
stocks that were coactive, such as the chart described herein with reference to Fig. 7. 

[0028] In yet another exemplary embodiment, a correlation map is plotted. To plot the 
correlation map, a correlation coefficient array is first generated for all of the stocks. The 
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correlation coefficients are defined as C(A,B)=number of times stock A and B are 
coactive divided by the number of times stock A is active. For this purpose, coactive 
means active at the same time, or within a specified number of time intervals of each 
other. The number of time intervals may be specified by a user or a default number such 
as one time increment may be used. The number of correlation coefficients will be equal 
to the square of the number of stocks observed. A correlation map is then drawn 
consisting of a map of all stocks with lines between each pair of stocks having a line 
thickness proportional to the correlation coefficient of those two stocks. An example of 
such a correlation map is illustrated in Fig. 7. There, an icon representing each observed 
stock 701, 703, 705, 707, 709, 711 is plotted around a circle 713. The thickness of line 
717 is proportional to the magnitude of the correlation coefficient for stocks 701 and 709. 
Line 715, which appears thicker than line 717, indicates that the correlation between 
stocks 705 and 709 is stronger than the correlation between stocks 701 and 709. 
Similarly, line 719, which appears thicker than lines 715 or 717, indicates that the 
correlation between stocks 701 and 705 is stronger than the correlation between stocks 
701 and 709 or stocks 705 and 709. If the correlation coefficient is below a 
predetermined threshold amount, the corresponding line may be omitted from the 
correlation map. The predetermined threshold amount may be specified by a user or a 
default threshold may be used. 

[0029] In still another exemplary embodiment, a cross correlogram is drawn to show 
potential causality among stock activity. This can be used to find stocks with events that 
consistently precede or follow events of another stock. A cross correlogram simply 
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creates a histogram of the time intervals between events in two specified stocks. A line of 
height proportional to the number of times the second stock is active one time interval 
following activity by the first stock is plotted at +1 on the x-axis of the histogram. A line 
of height proportional to the number of times the second stock is active two time intervals 
following activity by the first stock is plotted at +2 on the x-axis of the histogram, and so 
on. An example of such a cross correlogram is illustrated in Fig. 6. The line 601 
represents the number of occasions the first and second stocks were active at the same 
time, while line 607 represents the number of times the second stock was active three 
time intervals after the first stock was active. A cross correlogram may be plotted for a 
single stock to detect temporal characteristics in the stock's activity such as the fact that 
the stock is active with a period of every three time intervals a certain number of times 
during the period of observation. 

[0030] IDL code implementing all of the preceding steps of the exemplary embodiment 
is attached hereto as Appendix A. The procedure "MultiStock" and "MultiStock_event" 
are the main procedures. All relevant sub-procedures and functions are also included in 
Appendix A. 

[0031] An exemplary embodiment related to the cross-correlogram provides for 
displaying what is referred to as a "spike triggered average", which consists of a time 
series "movie" showing activity occurring in one or more stocks under investigation 
relative to activity in a selected stock. In this embodiment, a particular reference stock is 
selected. A data window consisting of a number of frames before and after events 
occurring in the selected stock (known as primary events) is then chosen or a default 
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number of frames may be used, such as ten. In the event ten frames are chosen, the 
resulting movie will consist of twenty-one frames, ten frames corresponding to the ten 
time periods before each event occurring in the reference stock, one frame corresponding 
to the time of each event in the reference stock and ten frames corresponding to the ten 
time periods after each event in the reference stock. 

[0032] Each frame of the movie will consist of a representation of all stocks under 
investigation. An example of such a frame is shown in Fig. 8. There, frame 800 consists 
of several icons 801, 803, 805, 807, 809 and 811, each corresponding to a stock under 
investigation. Each icon may be a solid square. The representations may also include 
ticker symbols 802, 804, 806, 808, 810 and 812 to further identify the stocks under 
investigation. A parameter of the icon for each stock, such as the color of the icon, is 
varied in each frame of the movie. The parameter varies in each frame to correspond to 
the frequency that events occur in the stock under investigation (known as secondary 
events) at the corresponding number of time periods before or after an event occurs in the 
reference stock. 

[0033] For example, if the reference stock selected had respective events at times t=20 
and t=50 and a movie length of twenty-one frames was selected, corresponding to ten 
frames before and ten frames after each primary event (i.e. an event in the reference 
stock), the movie would appear as follows. The first frame would be derived based on 
events occurring in the stocks under investigation at time t=10 and t=40 (i.e. 10 time 
periods before the respective events in the reference stock). Thus, if the first stock under 
investigation had an event at time t=10 and t=40, the icon parameter for that stock that is 
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displayed in the first frame would correspond to an event always occurring ten frames 
before an event in the reference stock, for example the icon color may be red. If the stock 
under investigation instead had an event at time t=10, but not at time t=40, the icon 
parameter for that stock that is displayed in the first frame would correspond to an event 
occurring half the time ten frames before an event in the reference stock, for example the 
icon color may be orange. The process is repeated for each stock under investigation for 
each of the frames in the spike triggered average movie. The resultant movie will 
illustrate the frequency that events occur in the stocks under investigation at the 
corresponding number of time periods before or after events occurring in the reference 
stock. This information may be used to uncover possible causality in the temporal 
domain among the stocks by identifying stocks whose activity appears to trigger or be 
triggered by activity in other stocks. 

[0034] In a still further exemplary embodiment, the data is analyzed in step 105 of Fig. 
1 by finding a hidden Markov state sequence from the event array. This embodiment 
uses the principal of Hidden Markov modeling described in Rabiner, A Tutorial on 
Hidden Markov Models and Selected Applications in Speech Recognition , Proceedings 
of the IEEE, vol. 77 pp. 257-286 (1989), which is incorporated by reference herein. 
Essentially, a Markov model is a way of modeling a series of observations as functions of 
a series of Markov states. Each Markov state has an associated probability function 
which determines the likelihood of moving from that state directly to any other state. 
Moreover, there is an associated initial probability matrix which determines the 
likelihood the system will begin in any particular Markov state. In a hidden Markov 
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Model, the Markov states are not directly observable. Instead, each state has an 
associated probability of producing a particular observable event. A complete Markov 
model requires the specification of the number of Markov states (N); the number of 
producible observations per state (M); the state transition probability matrix (A), where 
each element ay of A is the probability of moving directly from state i to state j; the 
observation probability distribution matrix (B), where each element bj(k) of B is the 
probability of producing observation k while in state i; and the initial state distribution 
(P), where each element pi of P is the probability of beginning the Markov sequence in 
state i. 

[0035] In the exemplary embodiment, it is assumed that the number of times events 
occur in a stock within each Markov state follows the Poisson distribution. Thus, each 
stock in each state has an associated Poisson Lambda parameter, which can be 
understood in the exemplary embodiment to correspond to the rate at which events occur 
in the stock. The set of all of these Lambda parameters is then assumed to be the B 
matrix. Given the estimations of the Markov Model parameters, the method uses the 
Viterbi algorithm to find the single best state sequence, i.e. the sequence of Markov states 
that most likely occurred to generate the observed results. The number of Markov states 
N may be selected by the user, or a default number such as six states may be used. The 
Viterbi algorithm is described as follows: 

Initialization: 

l<i<N, (1) 
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vi(0=o. ( 2 ) 



Recursion: 



5/y;-mg^ M &)B # ] 6,(0,) 2<i<T 
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l<i<N 



(3) 



(4) 
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Termination; 

z>* = max[8 r (/jl (5) 

1<*<AT v /y 

q* T = org max^> T (6) 



Path (backtracking): 

q>vM«) t = T-l,T-2,-,l. (7) 

[0036] In the algorithm, 5 t (i) represents the highest probability along a single path 
through all possible Markov state sequences up to time t that accounts for the first t 
observations (O t ) and ends in state i. \\f is used to store the argument which maximizes 
5 t (i). Once a possible state sequence q t * is generated, the state sequence plot can be 
generated such as the one shown in Fig. 9. In that example, six states are shown, 
corresponding to horizontal lines 901, 903, 905, 907, 909, 911. Each point on the plot 
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represents the Markov state the model is in at the relevant time. For example, point 913 
represents the Markov model being in state 903 while point 915 represents the model 
being in state 907. Each different state represents differing behavior of the stocks. For 
example, one group of stocks may exhibit events of interest more frequently than the 
remaining stocks when the model is in the first state 901, while those same stocks may 
exhibit fewer or no events when the model is in the second state 903. Correspondingly, 
another group of stocks may exhibit more frequent events of interest while in the third 
state 905 than other stocks and fewer events of interest while in the fourth state 907. 

[0037] A cross-correlogram between stocks in a selected state can be plotted using the 
methodology previously described, where only event data corresponding to the time the 
model is in the selected state is used in generating the cross-correlogram. The state may 
be selected by the user or a default state such as the first state may be used. 

[0038] IDL code implementing the preceding embodiment involving the hidden 
Markov model is attached hereto as Appendix B. The procedure "hidden_markov" and 
"hidden_markov_event" are the main procedures. All relevant sub-procedures and 
functions are also included in Appendix B. 

[0039] In a yet further exemplary embodiment the data is analyzed by performing a 
singular valued decomposition (SVD) on the data in the input stock data array, such as 
that shown in Fig. 3. In this embodiment, it is not necessary to detect events or store 
events in an event array. A singular valued decomposition takes advantage of the fact 
that in some sets of data produced from N different sources, such as N different stocks, 
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some of the stocks will not be creating independent data. In other words, there may be 
degeneracy in the data, which allows the data set to be decomposed into a number of 
eigenmodes i.e., orthogonal eigenvectors, with the eigenvalue (or singular value) 
representing the weight of the eigenvector in the system. 

[0040] In a singular valued decomposition, the data set is reduced from N dimensions, 
where N is the number of selected stocks, to d dimensions, where d is the specified 
number of eigenmodes and is less than N. The SVD algorithm, which is well known to 
one of ordinary skill in the art and is specified in the code in Appendix C, fits the 
observed stock data to a data model that is a linear combination of d number of functions 
of the spaces of data (such as time and stock price). Since d is specified rather than 
calculated by looking for degeneracy in the data, the resultant decomposition constitutes 
an approximation. Minimizing the sum of the squares of the errors in the approximation 
to the model, the SVD algorithm discards the eigenmodes corresponding to the smallest 
N-<2 eigenvalues. 

[0041] The stock data may be preprocessed before the SVD is performed by 
subtracting the median from each stock's closing price data. In other words, for each 
stock, a median is calculated and subtracted from each closing price entry for that stock. 
Additionally, when a positivity constraint is employed in the SVD algorithm (i.e. when 
only stock prices rising above the baseline are considered) an absolute value of the 
resultant data may be taken to ensure that downward events (i.e. drops in stock prices 
below the baseline) are considered in performing the SVD. 
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[0042] In this embodiment, the result that is plotted for visual analysis may be the level 
of each stock's contribution to each of the calculated d eigenmodes. For example, the 
result may be displayed in the format shown in Fig. 8, with each stock represented by an 
icon 801, 803, 805, 807, 809 and 811 and optionally a ticker symbol 802, 804, 806, 808, 
810 and 812. A parameter of the icon, such as its color, may be adjusted to represent the 
level of the stock's contribution to the displayed eigenmode. A separate plot can be 
generated for each of the calculated d eigenmodes. 

[0043] Alternatively, a plot, such as that shown in Fig. 10 may be generated to display 
the results of the SVD. This plot 1000, which displays singular values on the y-axis and 
mode number on the x-axis, represents the power of each mode in explaining the variance 
of the data set (i.e. the strength with which each of the calculated modes explains the 
tendency of the stock prices to deviate from the baseline). The example plot 1000 shows 
that most of the variance is explained by mode 0 (1006), mode 1 (1007) and mode 2 
(1008), while modes 3 (1009), 4 (1010) and 5 (1011) explain little of the activity in the 
data set. 

[0044] A third visualization useful to show the result of the SVD is shown in Fig. 11. 
In that example, three windows 1101, 1003 and 1005 are shown. The user first selects 
the mode for which data should be displayed, such as by using the slider bar 1 1 19. In the 
top window 1101, an icon for each stock (e.g. 1107, 1009) in the data set is displayed, 
with the stock's position on the y-axis corresponding to the strength with which that stock 
participates in the selected mode. The middle window 1103 shows a time series 
representation of the selected mode. In other words, window 1 103 displays the aggregate 
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stock activity corresponding to the selected mode. The bottom window 1105 is a 
superimposed plot of all of the stocks participating in the selected mode. As can be seen, 
the spike occurring around time day 300 (1115) in the bottom plot 1105 corresponds to 
the spike occurring at the same time (1111) in the aggregate mode activity shown in the 
middle plot 1103. Similarly, the spike occurring around day 480 (1117) in the bottom 
plot 1105 corresponds to the spike occurring at the same time (1113) in the middle plot 
1 103. Thus, it can be seen that activity in the identified stocks shown in the bottom plot 
1 105 does constitute the activity of the mode shown in the middle plot 1 103. 

[0045] DDL code implementing the preceding embodiment involving the singular value 
decomposition algorithm is attached hereto as Appendix C. The procedure "ssvd_gui" 
and "ssvd_gui__event" are the main procedures. All relevant sub-procedures and 
functions are also included in Appendix C. 

[0046] Although the present invention has been described in detail with reference to 
exemplary embodiments thereof, it should be understood that various changes, 
substitutions and alterations can be made hereto without departing from the scope or 
spirit of the invention as defined by the appended claims. 
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APPENDIX A 



pro choose_correl 

common mother__com, pixel_array, yes_no_yalues , coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common things_com, state3 

base4=WIDGET_BASE (/ROW, title= 1 Choose pvalue & windows i ze 1 ) 
fieldl=CW_FIELD (base4, / RE TURN_E VENTS , /FLOATING, TITLE='p value for chi 
squared 1 , $ 

VALUE = . 05 , UVALUE=0 ) 
f ield2=CW_FIELD (base4, / RETURN_E VENTS , / INTEGER, TITLE= 1 window size, ONLY 

ODD ' , $ 

VALUE=1, UVALUE=0) 
buttonl=WIDGET_BUTTON(base4, VALUE =' Go Cross Correlate', UVALUE=2) 
WIDGET_CONTROL, /realize, base4 

state3={f ieldl:f ieldl, f ield2 :f ield2, button :buttonl} 
WIDGET_CONTROL, WIDGET^INFO (base4 , /CHILD), S ET_UVALUE = s t a t e 3 

xmanager, ' choose_correl ' , base4 

end 

! NAME: — - — - " — 

; choose__correl_event , event 

; SYNOPSIS: 
,( j ; choose_correl_event , event 

CO ; DESCRIPTION: 

|ij ; This handles the events for choose_correl . Together they allow the user 

iii to 

|.| ; ; find the cross correlation coefficients and specify a p value and bin size. 

; As a result, the cross correlation matrix is printed to a file in the working 

directory 
J,;*' ; EVENT HANDLER 



B 



pro choos e_cor r e l_event , event 



; global variables_ 



common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
frame_no, f ile_name, $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common things_com, state3 

common flags, cells_def ined, spikes_def ined, single_plot_def ined, 
correl_coef_def ined 

stateholder3=WIDGET_INFO (event . handler , /CHILD) 
WIDGET_CONTROL, stateholder3 , GET_UVALUE=state3 
WIDGET_CONTROL, state3 . f ieldl , GET_VALUE=set_j>_yalue 
WIDGET_CONTROL, state3 . f ield2 , GET_VALUE=wins i ze 

; calculate the correlation coefficients for all the cell pairs 



if ( (set_p_yalue It 0) OR (set_p_value gt 1) ) then begin 

mess=WIDGET_MESSAGE ( 1 The p must be between 0 and 1!', /ERROR) 
endif else begin 

if (winsize ge frame_no) then begin 
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mess=WIDGET_MESSAGE ( 1 The window size cannot be larger than the total 
number of frames in the movie !' , /ERROR) 
endif else begin 

if ( (winsize mod 2) eq 0) then begin 

mess=WIDGET_MESSAGE (' Please enter an odd window size for 
uniformity', / INFORMATION) 
endif else begin 

coef =correl_coef f (yes_no_values , p=set_p_value, cells=cell_no, $ 

win=wins i ze , f rame_no ) 
close, 1 

openr, 1, 'correl.dat' 
stat_stuf f =f stat (1) 
f ile_size=stat_stuf f . size 
close, 1 

print, file_size 

/if (file_size gt 8112) then begin 

mess=WIDGET_MESSAGE ( 1 File is too large to display through a widget. 
Open it manually. If you did not save it, it is named correl.dat 1 , 
/ INFORMATION) 

; correl_coef_def ined=l 

; endif else begin 

xdisplayf ile, ' correl . dat 1 , title = "Correlation Coefficient Matrix", 
}rj group = event, top, width = 85, height = 45 
correl_coef_def ined=l 

f ile2=pickf ile (/write, f ile= ' Correl at ion_Coef f_exp_#_' ) 
}f f if (file2 eq 1 ' ) then begin 

jffj mess=WIDGET_MESSAGE ( 'The Correlation Coefficient Matrix has not been 

H saved. 1 , / INFORMATION) 
fjlj endif else begin 

CO get_lun, lunl 

* openw, lunl, file2 

m printf , lunl , coef 



e 



m 



free_lun, lunl 
close, lunl 
; endelse 
endelse 
endelse 
endelse 
endelse 

end 



pro choose_crosscorr 
base=Widget_base (/column) 

buttonl=widget_button (base, value= 1 Standard Cross Correlogram 1 , uvalue=l) 

button2=widget_button (base, value='Cell firing rate Cross 



Correlogram' ,uvalue=2) 
widget_control , base, /realize 



xmanager, ' choose_crosscorr ' , base 
end 

pro choose_cross cor r_e vent , ev 
widget_control , ev . id , get_uvalue=uval 
case uval of 

1 : draw_cross 

2 : draw_cross2 
endcase 
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end 



; NAME: 

choose_threshold 
; DESCRIPTION: 

; This program creates the widget with which the user can change the 

thresholds for the 

; detection of spikes in cells. It is called by calling Find Spikes to 
Multicell . 

; It uses MAKE__B I NARY (see choose_threshold_event ) to make the binary array, 
pro choose_threshold 

THE GLOBAL VARIABLES 



common mother__com, pixel_array , yes_no_values , coef , location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 

common with_choose_cells_com, thresholdl, threshold2 7 rms_threshold 

common things_com, state3 

base5=WIDGET_BASE (/Column, title= ' Choose threshold for spike detection') 
^ fieldl=CW_FIELD (baseB, / RE TURN_E VENTS , /FLOATING, TITLE= ' Threshold for 2 

consec . frames to be a spike ' , $ 
1% VALUE=threshl, UVALUE=0) 

J!;:: f ield2=CW_FIELD (base5, /RETURN_EVENTS, /FLOATING, TITLE= » Threshold for 3 
2;:^ consec> frames to be a spike 1 , $ 
W VALUE=thresh2 , UVALUE=0) 

H buttonl=WIDGET_BUTT0N(base5, VALUE- 1 Find Spikes ' , UVALUE=2) 
CO WIDGET_CONTROL, /realize, bases 

£8 state3={fieldl:fieldl, f ield2 : f ield2 , button : but tonl} 

3 WIDGET_CONTROL , WIDGET__INF0 (base5 , /CHILD), SET_UVALUE = s t a t e 3 

xmanager, 1 choose_threshold' , bases 

w. END 

m — 

; name: 

w ; choose_threshold_event 
fr A ; DESCRIPTION: 

/ This is the event handler for the choose_threshold widget. 

; With this, the user can change the threshold value used in determining which 

increases in 

; Calcium intensity correspond to action potentials.lt uses make_binary .pro to 
convert the 

; pixel or deltaF arrays into 0 or Is . 

EVENT HANDLER 



pro choose_threshold_event, event 

common mother_com, pixel_array, yes_no_values , coef, location, cell__no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common things_com, state3 

common with_choose_cells_com, thresholdl, threshold2 , rms_threshold 
common flags, cells_def ined, spikes_def ined, single_jplot_def ined, 
correl_coef_def ined 

stateholder3=WIDGET_INF0 (event . handler , /CHILD) 
WIDGET_C0NTR0L, stateholder3 , GET_UVALUE = s t a t e 3 
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WIDGET_CONTROL, state3 . f ieldl , GET_VALUE=threshl 
WIDGET_CONTROL , state3 . f ield2 , GE T__VALUE = t hr e s h2 

if ( (threshl eq 0) OR (thresh2 eq 0)) then begin 

mess=WIDGET__MESSAGE ( 1 Threshold cannot be OP, /Error) 
endif else begin 

WIDGET_CONTROL, event . top, /hourglass 
for cell = 0, cell_no - 1 do begin 
thresholdl [cell] = threshl 
threshold2 [cell] = thresh2 

endf or 

yes_no_values=make_binary (pixel_array , f rame_no, cell_no, 
thresholdl , threshold2 ) 

spikes jief ined=l /this lets other programs 

now that spikes have been defined 

mess=WIDGET_MESSAGE ( 'Done finding spikes boss. 1 , / INFORMATION) 
endelse 

end 

contingency_values 
DESCRIPTION: 

This program creates a contingency table for the data of two spike trains 
and from it 
" • gives a chi squared value. 



METHOD: observed_array is an array of 4 elements they are as follows: 

observed_array (0) ^number of hit_hits, i.e the number of times the first 



rjj 

\| cell and the second 

m 



fp; ; cell have spikes at exactly the same time. 



5 : 



observed_array (1) = number of observed hitjnisses 
observed_array (2) = number of observed miss_hits 
observed_array (3) = number of observed miss_misses 



This is called by CORREL_COEFF . PRO to calculate the correlation coefficients 
j*' between all the cells of 
f!j ; an analysis. Using the contingency table, we get a chi squared value which we 
O compare with the chi squared 

14; ; value created when the user specifies the p for the data to be significant. 
Each pair of cells has its 
; own contingency value. 

function contingency_values , filearrayl, filearray2 

;this program will print out a contingency table and then print out the chi- 
square value. 

;it also print out the probability that the null hypothesis is true 
;Null hypothesis: The two files are independant 

; Constraint: THIS WILL ONLY WORK IF THE SPIKE TRAINS FOR THE CELLS ARE IN TERMS 
OF 0 AND 1 

/declaring the arrays 

observed_array=f Itarr (4) 

expected__array=f ltarr (4) 

result=f ltarr (4) 

/This makes a spike for celll = 2 where a spike for cell2=l. Therefore 
celll-cell2 will be 0 only 
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;at times when both celll and cell2 have a 0 (at miss_miss) . A 
;Due to this, a hit_miss can be found by seing where the above difference 
array is equal to 1 . 

/eg. , 100101 -> 200202 

filearrayl=2*f ilearrayl 

;the following array is used to calculate all of the values needed in the 

table 

;we subtract the spike trains of the two cells into dif f erence_array . 
/therefore if dif f erence_array has a 0, there 

dif ference_array=f ilearrayl -filearray2 
/this is the m(l,l) values or the hit_hit 

;a spike in f ilearrayl will be a 2 . That in filearray2 will be a 1. 
Therefore, in the difference array 

; a hit_hit will be represented by a (2-1) or a 1 
arr_hit_hit=dif f erence_array eq 1 
if (total {arr__hit_hit) gt 0) then begin 

observed_array (0) =total (arr_hit_hit) 
f endif else begin 

P observed_array ( 0 ) = 0 

I'll endelse 
\ a 



P 



;this is the m(0,0) values or miss_miss 
arr_temp3=dif f erence_array eq 0 

if (total (arr_temp3) eq 0) then begin 

observed array (3) =0 
endif else begin 

miss_miss_data=total (arr_temp3) 
Hj observed_array (3) =miss_miss_data 

p endelse 



m 

y, 



0))) 



;observed_array (3) =n_elements (dif ference_array (where (dif f erence_array eq 



;this is the 111(1,0) values or the hit_miss, i.e the number of times when 
at an instant in time, 

;the first cell has a spike, but the second cell does not. 
arr_hit_miss=dif f erence_array eq 2 

if (total (arr_hit_miss) eq 0 ) then begin 

observed_array (1) =0 
endif else begin 

hi t_miss_data= total (arr_hit_miss) 
observed_array (1) =hit_miss_data 
endelse 

/this is the m(0,l) values or the miss_hit, i.e the number of times when 
at an instant in time, the 

/second cell has a spike while the first cell does not. (celll (t)=0, but 
cell2 (t) -1) 

arr_miss_hit=dif f erence_array It 0 

if (total (arr_miss_hit) eq 0) then begin 
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observed_array (2) =0 

endif else begin 

mi s s_hi t_da t a= t ot al ( arr_mi s sjii t ) 
observeci_array (2) =miss_hit__data 

endelse 



sum_f irst_row=observed_array (0) +observed_array (1) 
sum_second_row=observed_array (2) +observed_array (3 ) 
sum_f irst_col=observed_array ( 0 ) +observed_array (2 ) 
sum_second_col=observed_array ( 1 ) +observed_array ( 3 ) 

n__t o t a 1 = s um_f i r s t_r ow+ s um_s e c ond__r ow 

; now we calculate the observed and expected values for chi squared 
;for j=l, i=l 

expected_array (0) = (sum_f irst_row*sum_f irst_col) /n_total 
;for i=l, j=0 

expected_array (1) = (sum_f irst_row*sum_second_col) /n_total 
;for i = 0, j =1 

expected__array (2) = (sum_second_row*sum_f irst_col) /n_total 
;for i=0, j=0 

expected_array (3) = (sum_second_row*sum_second_col) /n_total 



for k=0 / 3 do begin 

result (k) = ( (observed_array (k) -expected_array (k) ) * (observed_array (k) - 
FU expected_array (k) ) ) / expect ed_array (k) 
%i endfor 

answer= total ( result ) 
fp return, answer 

7 end 

|:f#: — — — 

! a ; NAME: correl_coef f . pro 

pi ; INTRODUCTION: 

H ; This routine is used to output a correlation coeff . matrix with the 

«[:f insignificant 

S3 ; correlation coeff =0. It uses the functions Contingency__values and 

|4 f ind_jnatches 

INPUTS : 

p_value-> This is the probability that a random variable is greater than 
the cut off value in 

; a chi squared distribution. For example, p=.05 means the cut off value should 
be high enough 

so there is only a 5 percent chance that any random number can be greater. 

cell__number-> the number of cells chosen from the slice. 

winsize-> the size of the window to look for coincident spike in. The 
default is 1 



DESCRIPTION: 

The function does the following: 

1. Uses the IDL function chi_sqrcvf to find the cut off value for 
the input p value. It is 

; called as chi_sqrcvf (p_value , degrees of freedom) 
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2. Creates a contingency table for a pair of cells by the home-made 

function 

; contingency_values . pro. From this table, it then finds the chi squared value 
for the null hypothesis that 
; the cells are independant. 

3. It prints out all of the correlation coefficients that are 
significant (fail the null 

; hypothesis, but pass the cut_off test) and indicate the cells of the pair are 
dependant . 

; OUPUT: 

a matrix of the filtered correlation coefficients. 
; NOTE : To see all of the cross correlation coefficients, make the p_value=ll 

function correl_coef f , yes_no_values , p=p_value, cells=cell_number , win=winsize, 
f rame_no 

forward_f unction contingency_values , find_matches 
output_array=f ltarr (cell_number, cell_number) 



-FIND THE CUT OFF VALUE - 



limit=chisqr_cvf (p_value, 1) ; the degree of freedom is 1 in our 

H case 

ru 

; t ; CONTINGENCY TEST AND CORREL COEFF CALC. 

H _ 

Ul bin_no=f rame_no/winsize ;using the window size, this is the 

* number of bins 

M- lef trover s=frame_no mod winsize /these are the frames that 

}^ didn't fit in the bins, but are still in the array (at the end) 
in total elements=bin no+left overs ;this is the total number of 

ft; elements in the binned array 

if (total__elements It 1) then begin 

mess=WIDGETJVIESSAGE ( 'Window size too large or two few frames.', /ERROR) 
r ' endif else begin 

close, 1 

IF (winsize eq 1) THEN BEGIN ; IF THE WINDOW SIZE IS 

ONE, NO BINNING IS DONE! 

for i=0, cell_number-l do begin 

for j=0, cell_number-l do begin 

temp_value==contingency_values (yes_no_values (i, *) , 

yes_no_values ( j , *)) 

if (temp_yalue ge limit) then begin 
output_array ( i , 
j) =f ind_matches (yes_no_values (i, *) , yes_no_values ( j , *) , $ 

winsize) 

endif 

endf or 

endf or 

ENDIF ELSE BEGIN 
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cell l_binned= intarr ( tot al_element s ) 
cell2_binned=intarr (total_elements) 
stop_index=-l 

index=0 ^ n case 

no bins are made, then the stop index used in the second for loop below will be 
used 

for i=0, cell_number-l do begin 

for j=0, cell_number-l do begin 

for bin_number=0, bin_no-l do begin 
; making the binned cell array 

start_index=bin_number*winsize 
;this is the starting point for a bin 

s t op_index= s t art_index+wins i ze - 1 
;this is the ending point for a bin 

celll_binned (bin_number) =total (yes_no_yalues (i, start_index: stop_index) ) 
;we are binning the input arrays for the contingency table 

cell2_binned (bin_number) =total (yes_no_values ( j , start_index : stop_index) ) 
; f or the second cell being looked at 
endf or 
index=0 

H for bin_number=stop_index+l, frame_no-l do begin 

O ;add the elements that don't fit in the bins 

0 celll_binned (bin_no+index) =yes_no_values (i , bin_number) 

p cell2 jDinned (bin_no+index) =yes_no_values ( j , bin_number) 

| ; y i ndex= i ndex+ 1 

endf or 

celll_binned=celll_binned ge 1 
;this converts the arrays into binary ones (eg if three elements in a 
™ bin each have a value 1 then the total will be 3 that will be converted to 1 
* here) 

H' cell2_binned=cell2_binned ge 1 

M temp_value=contingency_values (cellljainned, cell2 _binned) 

0 if <temp_value ge limit) then begin 

jfjj output_array ( i , j ) =f ind_matches (yes_no_values (i, 

fi *) yes no values ( j , *) , $ 



winsize) 
endif 

endfor 

endf or 
ENDELSE 
close, 1 

openw, 1, ■ correl.dat' ; CHOOSE_CORREL USES XDISPLAYFILE AND THIS 

FILE TO WRITE 

printf , 1 , output_array 

close, 1 

return, output_array ;this is the correl coeff matrix 

endelse 

end 



NAME: 



CORREL MAP IMAGE PLANE 



PURPOSE : 

This procedure creates a circular representation of correlation between 
cells of the slice. 
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Those cells that 

are correlated in either direction are joined by lines which are 
proportional to their 
; correlation coefficients 

pro correl_map_image_jplane, coef_array 

common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
f rame_no, f ile_name, $ 

total_f rame^no, time_resolution, x_size, y_size, box_size 

for s=0, cell_no-l do begin 

; ; 1) draw boxes over active cells 

;; we only worry about the active cells- cells that have at least 

one spike 

if total (yes_no_values (s, *) ) gt 0 then begin 

x=location(s) .coord (0) ;the x axis of where you clicked on 

cell 

y=location(s) .coord (1) ;the y axis of where you clicked on 

cell 

for j=-box_size, box_size do begin 

plots, [x+j, x+j], [y-box_size, y+box_size] , /device 

endf or 

xyouts, x-50, y, s+1, font=-l, charsize=l . 9 , /device 



U ;; 2) draw lines between 'correlated' cells 

Q for t=0, cell_no-l do begin 

p if (coef_array (s, t) gt 0) then begin 

ffl plots, 



; y [location (s) . coord (0) , location (t) . coord (0) ] , [location (s) . coord (1) , 
m location (t) . coord (1) ] ,$ 

thick=2*coef_array (s, t) , /device 

^ endif 
* endfor 
H endif 
H endfor 
Q !p.font=0 
fU end 

C9 — 

u 

; Name: count_random_hits_2_manyX, least_no_of_matches 

; Synopsis : count_random_hits_2_manyX / random_array , least_no_of_matches 

; Description: This program tests to see how many times more than 

least_no_ofjnatch.es two random cells 

; fire simultaneously. 

; Description of variables : 

; yes_no_values : the binary data 

random_array: a random array created by RAND0M_TEST looking at the spikes 
in yes_no_values 

least_no_of _matches : Minimum number of times that two cells must fire to 
be considered a coactive pair 

; window_size: number of frames before and after to look for a coincident 

spike in the other cell 

function count_random_hits_2_manyX / random_ar ray , least_no_of_matches , 
window_s i z e , num_a c t i ve_c e lis 

common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
f r ame_no , f i 1 e_name , $ 
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total_f rame^no, time_re solution, x_size, y_size, box_size 
random_hits=0 

for cell_l_counter = 0, num_active_cells - 2 do begin 

for cell_2_counter = cell_l_counter + 1, num_active_cells - 1 do 

begin 

temp_hits = 0 

for window_counter = -window_size , window_size do begin 
temp_cell_l = intarr (f rame_no + (2 * window_size) ) 
temp_cell_2 = intarr (frame_no + (2 * window_size) ) 
temp_cell_l (window_size : f rame_no+window_size-l) = 

random_array (cell_l_counter, *) 

temp_cell_2 (window_size : f rame_no+window_size-l) = 

random_array (cell_2_counter, *) 

temp_cell = (temp_cell_JL * shif t (temp_cell_2 , 

window__counter) ) 

temp_hits = temp_hits + total (temp_cell) 

endf or 

if (tetnp_hits ge least_no_of_matches) then begin 
random_hits = random_hits + 1 

endif 

f ™ end for 

W end for 

D 

Q return, random_hits 

III end 

^| ; Name: count_random_matches 

13 ; Description: this program takes in a random array of spike trains and counts 

p the number of sets of coactive spikes 

" ?> ; (where the number of spikes in a set is ' times_repeated » , e.g. if 2, we 

y count the number of all possible pairs 

f" ; of coactive spikes 

iJ ; Defintion of variables: 

fU ; random_array is the array of the location of spikes for the simulated data 

Qj ; times_repeated is the number of times the spikes are clustered together in 

Mi calculating the 

combination. For example, if looking at pairs of two, 
times_repeated=2 . If looking for cells firing 3X together, it is 3. 

Formula Used: C(n,r)= n!/(r! * (n-r) ! ) 

function count_randomjnatches, random_array , times_repeated, window_size 

common mother__com, pixel_array, yes_no_values , coef , location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 

;find the number of hits for the simulated data. 

; we can do this by summing the rows (cells) of one column (frame) since 
spikes are 

; represented by 1 and then finding the combinations of this taken 
times_repeated at a time 

; eg. cell 1 : 001100110101010100101010 

cell 2: 001000100001000001000001 
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; If in frame 23, only cells 2, 3, and 8 fire, the sum of the frames will 
; be 3 and then the combination taken two at a time is C(3,2)=3 



random^matches = 0 

if (window_size eq 0) then begin 

for frames = 0, frame_no - 1 do begin 

matches = total ( random_array (* , frames)) /finds n or the 

number of hits at this frame 

if ((matches gt 0) and (matches ge times_repeated) ) then begin 
; ; find number of possible pairs 
matches _pair = factorial (matches) / 
(factorial (times_repeated) * factorial (matches -times_repeated) ) 
endif else begin 

matches_pair = 0 ; cant find factorial 

of a negative number 

endelse 

randomjnatches = random_matches + matches_pair 

endf or 
endif else begin 

; ; first compress spike trains according to window size 
binned_randon_array = intarr (n_elements (random_array [* , 0] ) , frame_no 
■;^ r - (window_size * 2)) 

S for f rame__counter = window_size, frame__no - (window__size + 1) do 

fcl begin 

p for window_counter = -window_size, window_size do begin 

fjj for cell_counter = 0, n_elements (random_array [* , 0] ) - 1 

\j do begin 

m if (random_array [cell_counter, f rame_counter + 



window_counter] ge 1) then begin 



b inned__r andom_ar r ay [eel l_count e r , 



! ; f rame_counter - window_size] = 1 
H 1 " endif 

|f«? endfor 
Q endfor 
Hi endfor 

pi ; ; then look for matches in the binned spike trains (although a 

| f | : single spike can fall into many bins) 

for f rame_counter = 0, n_elements (binned_random_array [0 , *] ) - 1 do 

begin 

matches = total (binned_random_array (* , f rame_counter) ) 
; finds number of hits at this frame 

if (matches ge times_repeated) then begin 
matches_pair = factorial (matches) / 
(factorial (times_repeated) * factorial (matches - times_repeated) ) ; finds 
number of possible pairs 

endif else begin 

matches_pair = 0 ;cant find factorial 

of a negative number 

endelse 

randomjnatches = randomjnatches + matches__pair /increases 
the matches count for the random data by that of this frame 
endfor 
endelse 

return , random_match.es 

end 
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; Name: count_random_matches_many 

; Description: this program takes in the actual binary spike train and the 
random one and sees how 

; the actual one compares to the random one. To test the significance, it looks 
at the number of 

; coincident spikes boths sets of data have and creates a uniform distribution 
for the random data 

; to see where the actual data lies. 

; Defintion of variables: yes_no_values is the binary spike train of the true 
data 

random_array is the array of the location of spikes for the 

simulated data 



function count_random_matches_many, random_array, times_repeated, ge_or_eq_test , 
window_size 

common mother__com, pixel_array, yes_no_values , coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_frame_no, time_resolution, x_size, y_size, box_size 

, , random matches = 0 

^ if (window_size eq 0) then begin 

W for frames = 0, frame_no - 1 do begin 

p matches = total (random_ar ray (* , frames)) ; finds number 

flj of hits at this frame 

; y if ( ( (ge__or_eq_test eq 0) and (matches ge times_repeated) ) or 

|f| ( (ge_or_eq_test eq 1) and (matches eq times_repeated) ) ) then begin 
J ? ! random_matches = random_matches + 1 

endif 

! ; , endf or 

j7 endif else begin 

W ;; first compress spike trains according to window size 

P binned_random_array = intarr (n_elements (random_array [* , 0] ) , frame_no 

flj - (window_size * 2)) 

f]j for f rame_counter = window_size, frame_no - (window_size + 1) do 

U begin 

for window_counter = -window^size, window__size do begin 

for cell_counter = 0, n_elements (random_array [* , 0] ) - 1 

do begin 

if (random_array [cell_counter, f rame_counter + 
window_counter] ge 1) then begin 

binned_random_array [eel l_count er , 

f rame_counter - window_size] = 1 

endif 

endf or 

endf or 

endf or 

; ; then look for matches in the binned spike trains (although a 
single spike can fall into many bins) 

for f rame_counter = 0, n_elements (binned_random_array [0 , *] ) - 1 do 

begin 

matches = total (binned_random_array (* , f rame^counter) ) 
/finds number of hits at this frame 

if ( ( (ge_or_eq_test eq 0) and (matches ge times_repeated) ) or 
( (ge_or_eq_test eq 1) and (matches eq times_repeated) ) ) then begin 
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random_matcb.es = random_matches + 1 

endif 

endf or 
endelse 

return, randomjnatches 

end 

Name: delete_spikes_widget 

Description: allows me to delete from all cells a spike in a particular frame. 

pro delete_spikes_widget 
common delete_spikes , statel2 

basel2=WIDGET_BASE (/Column, title=' Spike Deletion' ) 
buttonll=WIDGET_BUTT0N(basel2 / VALUE= ' Delete All Trailing Spikes', 
UVALUE=11) 

fieldl2=CW_FIELD(basel2, / RE TURN_E VENTS , / INTEGER , TITLE= T Frame for Global 
Spike Deletion', VALUE=0, UVALUE=0) 

buttonl2=WIDGET_BUTT0N(basel2, VALUE= ' Delete Spike in This Frame', 
UVALUE=2) 

fieldl3 = CW_FIELD(basel2, / RETURNEE VENT S , /FLOATING, TITLE= 1 Threshold in 
dF/F Units for Global Spike Deletion', VALUE=0.0, UVALUE==0) 



buttonl3=WIDGET BUTTON (base 12 , VALUE= ' Delete Spikes Less Than Threshold', 



□ 

j:f UVALUE=2) 

P WIDGET_CONTROL, /realize, basel2 

FU statel2={buttonll :buttonll, f ieldl2 : f ieldl2 , buttonl2 :buttonl2 / 

H| fieldl3 :fieldl3, buttonl3 :buttonl3 } 

(Q WIDGET_C0NTR0L, WIDGET_INF0 (basel2 , /CHILD), SET_UVALUE=statel2 

fp. 

** xmanager, ' delete_spikes_widget 1 , basel2 

j . end 

|;;|; 

m 
f f % 



Name : delete_spikes_widget_event 

Description: allows me to delete from all cells a spike in a particular frame. 



pro delete_spikes_widget_event , event 

common mother_com, pixel_array, yes_no_yalues, coef , location, cell_no, 
f r ame__no , f i 1 e__name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common delete_spikes , statel2 

stateholderl2=WIDGET_INF0 (event . handler, /CHILD) 
WIDGET_CONTROL, stateholderl2 , GET_UVALUE= S t a t e 1 2 
WIDGET_CONTROL, statel2 . f ieldl2 , GET_VALUE=f rame_to_kill 
WIDGET_CONTROL, statel2 . f ieldl3 , GET_VALUE=spike_unit__threshold 

; ; deleting ALL trailing spikes in ALL cells 
if (event. id eq statel2 . buttonll) then begin 

for cell_counter = 0, (cell_no - 1) do begin 

for f rame_counter = (frame_no - 1) , 2, -1 do begin 

if ( (yes_no_values [cell_counter , f rame_counter] eq 1) 
and ( (yes_no_yalues [cell_counter , frame_counter - 1] eq 1) or 
(yes_no_values [cell_counter , f rame_counter - 2] eq 1) ) ) then begin 
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yes_no_values [cell_counter , f rame_counter] = 0 

endif 

endfor 

endfor 

endif 

; ; delete spike in every cell that is in a given frame 
if (event. id eq statel2 -buttonl2) then begin 

if (frame_to_kill ge 0) and (f rame_to_kill It frame_no) then begin 
for cell_counter = 0, (celljio - 1) do begin 

yes_no_yalues [cell_counter , f rame_to_kill] = 0 

endfor 

endif 

endif 

; ; delete all spikes less than a raw dF/F threshold of intensity 
if (event. id eq statel2 .buttonl3 ) then begin 

for cell_counter = 0, (cell_no - 1) do begin 

for f rame_counter = 0, (frame_no - 1) do begin 

if (pixel_array [cell_counter , f rame_counter] It 
spike_unit_threshold) then begin 

yes no values [cell counter, frame counter] = 0 

|:f4; — — — — 

endif 

*jy endfor 
P endfor 
0 endif 

rtj end _ __ 

% '4 ; Name: draw_3D_plot 

m ; 

|p ; Description: draws all frames of all cells, with intensity of dF/F on the z- 
axis; 

a . ; i.e. it makes a pretty cool plot! 

}* ; 

™? pro draw_3D_plot , pixel_array 

!1J 

0 window, 4, title='All Stocks -3D 1 

M: surface, transpose (pixel_array) , /horizontal 

end 

; NAME: draw_cgram 

; 

; INTRODUCTION: 

This procedure is used to draw cross correlograms between any two cells. A 
cross correlogram 

; is a histogram of the time intervals between the two cells. If there is a peak 
at 0, the cells are 

; very highly correlated. This draws a cross correlogram between cell 2 with 
respect to cell 1. That is 

; a peak at +1 means that cell 2 leads cell 1 by 1 time frame. 
; 

; INPUTS : 

; pixel_array. This is the array created in either pixel_binary_cells or 

pixel_binary_box 

; depending on the choosing method you used. It contains the pixel values of the 
cells and their 

; converted binary form. It is the concated array of: 
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pix__array=intarr (cell_no, frame_no) {actually called 
pixel__array in those pro} 

. yes_no_values=intarr (cell__no, frame_no) 

where cell_no is the number of cells you chose 

and frame_no is the number of frames you chose to analyze 

therefore, what we are talking about is: 

pixel_array =intarr (cell_no+cell+_no, frame_no) 

to find the data on cell x, 
; pixel_array (x-1, *) will give the pixel values over time (from 

pix_values) 

pixel_array (cell_no+x-l, *) gives 1 for a peak else 0 (from 

ye s_no_va lue s ) 

; cell__numberl: This is the number of the first cell you want to cross 

correlate from 

cell_number2 : This is the number of the second cell you are considering 
■ . ; binsize: This is the window size you are calculating with. If any spikes fall 
f) * in this, they are 

Jjjf ; considered coincident. Default is 1 

y ; cell_no: Total number of cells you have chosen 

0 ; time: The time between two frames shot by the CCD camera (input from 

111 widget) 

" : f..fj * 

Cp ; DESCRIPTION: 

; This program looks at the yes_no_values part of pixel_array (look above) 

and converts the 

?, ; binary array into an array holding the positions of the spikes. It uses the in 
f7 built function 

; where () to do this. Then you use the function FIND_D I STANCE to calculate the 
M time differences between 

llj ; the spikes. It's source code is located in the same directory as this, 'All 
13 you need' . Briefly what 

H ; it does is it subtracts the two time arrays 

pro DRAW_CGRAM, yes_no_values , cell_numberl , cell_number2 , bin_size, cell_no, 
time, xmin, xmax 

/this procedure will make a correlogram between the two files input. 
;the functions it uses are : findgen (size) *bin_size+mintemp 
;openf ile, f ind_distance 

Ip .multi=0 

forward_function F IND__D I S TANCE 

write=string(format=' ( "cross -correlogram of cell",/, (iO) ) ', cell_numberl) 

write2=string (formats ' (" vs cell ",/,(i0))\ cell_number2 ) 

timel=where (yes_no_values (cell_numberl-l, *) gt 0) ;this finds where 

the spikes are 

time2=where (yes_no_values (cell_number2-l, *) gt 0) 
IF ((timel(O) eq -1) OR (time2 (0) eq -1)) then begin 

mess=WIDGET__MESSAGE ( ' One of the cells has no spikes') 
ENDIF ELSE BEGIN 
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dist emp = F IND_D I S TANCE (timel, time2) 

mintemp=min (distemp) 

maxtemp=max (distemp) 

xmin=f loat (xmin) /time 

xmax=f loat (xmax) /time 

size=xmax-xmin 

IF bin_size It 1 THEN BEGIN 

mess=WIDGET_MESSAGE ( 1 Binsize must be >=1 * ) 
END IF ELSE BEGIN 

hist=histogram (distemp, binsize=bin_size / min=xmin / max=xmax) 
plot, (findgen(size) *bin_size+xmin) *time / hist, psym=10, 
title=write (0) +write (1) +write2 (0) +write2 (1) , $ 

xtitle= ? time in seconds. 1 , ytitle='no of spikes' 

ENDELSE 
ENDELSE 

end 



pro DRAW_CGRAM2 , yes_no_values , numberl, number 2 , binsize, cell_no, time, xmin, 
xmax 

;this procedure will make a correlogram between the two files input. 
9 ;the functions it uses are : findgen (size) *bin_size+mintemp 
O ;openfile, find_di stance 

o 

fjj !p.multi = 0 

\J 

write=string(format=' ( "cross -correlogram of cell " , / , ( iO ) ) * , numberl) 
write2=string (formats (" vs cell ",/,(i0))\ number2) 

timel=where (yes_no_values (numberl -1, *) gt 0) /this finds where the 
spikes are 

time2=where (yes_no_values (number2 -1, *) gt 0) 
M' IF ((timel(O) eq -1) OR (time2 (0) eq -1)) then begin 

Q mess=WIDGET_MESSAGE ( ' One of the cells has no spikes 1 ) 

fjj END IF ELSE BEGIN 



U N_l=total (yes_no_values , 2 ) 

N=N_1 [numberl - 1] ;this finds the total number of spikes for 

cell numberl 

s i z eo f yn_va 1 = s i z e ( y e s__no_value s ) 
actual_f rame_no=sizeofyn_val [2] 

Si=intarr (actual_f rame_no) 
for t=0, actual_f rame_no - 1 do begin 
Si [t] =yes_no_values [numberl -1, t] 
endf or 

Sj=intarr (actual_frame__no) 
for t=0, actual_f rame_no - 1 do begin 
Sj [t] =yes_no_values [number2 -1, t] 
endf or 

nlim=fix(- actual_f rame_no / binsize) ;+l ;be careful - some data is lost 
when tha actual fr.# is not 

plim=f ix(actual_frame_no / binsize) ; -1 ;be careful -divisible by 
binsize 
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; print , 



'nlim - lower bound = 



nlim, 



plim - upper bound = ', plim 



R=fltarr{2, plim - nlim +1) 
index=0 

for l=nlim / plim do begin 
z = 0 

for t=0 / actual__f rame_no -1 do begin 
x=0 

for k=t+l*binsize, t+ (1+1) *binsize -1 do begin ;be careful with the -1 

if (k ge 0) and (k It actual_f rame_no) then begin 

x=x+Sj [k] 

endif 

endf or 

if {x gt 0) then begin 

x=l 

endif 

z=z + Si[t]*x ; (x is either 0 or 1 depending whether cell j fired in a 



xminl=f ix (xmin/ time) 
xmaxl=f ix (xmax/time) 

IF binsize It 1 THEN BEGIN 

mess=WIDGET__MESSAGE ( 'Binsize must be >=1') 
ENDIF ELSE BEGIN 

IF (center + xminl It 0) or (center + xmaxl gt index -1) then begin 

mess=Widget_Message ( 'X is out of range 1 ) 
Endif else begin 

xaxis- [R[0, center + xminl] * time , (INDGEN ( - xminl + xmaxl) 
+R [0, center+xminl] +1) *time] 

yaxis=fltarr (xmaxl -xminl +1) 

for i=center + xminl, center + xmaxl do begin 

yaxis [i - center - xminl] =R [1 , i] 

endfor 




given time bin) 
endfor 



R[0, index] =1 
R [1, index] =z 
index=index+l 
endfor 



for i^O, index -2 do begin 
R[l,i]=R[l,i]/(binsize * N) 
/print, R[l,i] 
endfor 



center=-l 

for i=0, index-2 do begin 

if (R[0,i] eq 0) then center=i 

endfor 
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plot, xaxis, yaxis , psym=10, 
title-write (0) +write (1) +write2 (0) +write2 (1) , $ 

xtitle=' time in seconds.', ytitle= ' firing rate 1 



endelse 
ENDELSE 
ENDELSE 

end 



if: 

m 

U 



pro draw_cross 

common mother_com, pixel__array, yes_no__values , coef , location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_frame_no, time_resolution, x_size, y_size, box_size 
common markov2 , yes_no_values3 , yes_no__temp /added by pvc3 

base=WIDGET_BASE (/column, title=' Cross Correlogram' ) 

draw=WIDGET_DRAW (base , xsize=400, ysize=400) 

number 1=CW_FIELD (base # title= f Enter Cell Number 1', VALUE=1, UVALUE=2,$ 

/ RETURNEE VENTS , / INTEGER ) 
number2=CW_FIELD (base, title-' Enter Cell Number 2', VALUE=2 , UVALUE=3,$ 

/ RE TURN_E VENT S , / INTEGER ) 
binsize=CWJFIELD(base, title=' Enter the binsize', VALUE=1, UVALUE=4,$ 

/RETURN_EVENTS , / INTEGER) 
xmin=CW__FIELD(base, title='xl=Lower bound', VALUE=-200, / RETURN_E VENTS , 
/INTEGER) 

xmax=CW_FIELD(base, title= 1 x2=Upper bound 1 , VALUE-200, / RE TURN_E VENTS , 
/INTEGER) 

text=WIDGET_TEXT (base, VALUE=string (' Total number of cells you chose:', 
cell_no) ) 

WIDGET_C0NTR0L , base, /realize 

holder = {number 1 : numberl , number2 : number2 , binsize : binsize , draw : draw, 
xmax : xmax , xmi n : xmi n } 

WIDGET_CONTROL, WIDGET_INFO (base , /Child), SET_UVALUE= holder 
xmanager, ' draw_cross 1 , base 
end 



pro draw_cross_event , event 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 

common markov2 , yes__no_values3 ; added by pvc3 

stateholder5=WIDGET_INF0 (event . handler , /Child) 
WIDGET_CONTROL, stateholder5 , GET_UVALUE = s t a t e 
WIDGET_CONTROL , state . draw, GE T_VALUE = w i n_i d 2 
WIDGET_CONTROL, state . numberl , GE T_VALUE = numbe r_ 1 
WIDGET_CONTROL , state . number2 , GET__VALUE = numb e r__2 
WIDGET_CONTROL, state .binsize , GET_VALUE=bin_size 
WIDGET_CONTROL , state . xmax , GET_VALUE=xmax 
WIDGET_CONTROL , state . xmin, GET_VALUE=xmin 



IF ( (number_l gt cell_no) OR (number_l It 1) ) OR ( (number_2 gt cell_no) OR 
(number_2 It 1) ) THEN BEGIN 

mess=WIDGET_MESSAGE (string (' Invalid cell number entered. You only 
chose 1 , byte (celljrio) , ' cells'), /ERROR) 
END IF ELSE BEGIN 
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wset, win_id2 /changed by pvc3 

draw_cgram, yes_no_values3 , numberJL, number_2 , bin_size, cell_no, 
time_re solution, $ 

xmin , xmax 

ENDELSE 

end 



pro draw_cross2 

common mother_com, pixel_array, yes_no_values, coef , location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box__size 
common markov2, yes_no_values3 , yes_no_temp /added by pvc3 

base=WIDGET_JBASE (/ column, titled Cell firing rate Cross Correlogram' ) 

draw=WIDGET_DRAW (base, xsize=400, ysize=40 0) 

numberl=CW_FIELD(base,title= 'Enter Cell Number l 1 , VALUE=1, UVALUE=2 , $ 

/ RE TURN_E VENTS , / INTEGER) 
number2==CW_FIELD(base, title= f Enter Cell Number 2', VALUE=2 , UVALUE=3,$ 

/ RE TURN_E VENTS , / INTEGER) 
binsize=CWJFIELD(base, title= 'Enter the binsize', VALUE=1, UVALUE=4,$ 
M : / RE TURN__E VENTS , /INTEGER) 

0 xmin=CW_FIELD(base, title- » xl=Lower bound 1 , VALUE=-200, / RETURN_E VENTS , 

p /INTEGER) 

5 xmax=CW_FIELD(base, title= 'x2=Upper bound 1 , VALUE=200, / RETURNEE VENTS , 

jf|"s / INTEGER) 

ill text=WIDGETJTEXT (base, VALUE=s t ring ( 1 Total number of cells you chose:', 

cell_no) ) 

^ WIDGET_CONTROL, base, /realize 

tt holder={numberl:numberl, number 2 : number 2 , binsize :binsize, drawrdraw, 

5 xmax: xmax, xmin: xmin} 

M : WIDGET_CONTROL, WIDGET_INFO (base , /Child), S ET_UVALUE =ho 1 de r 

M ; xmanager , 1 draw_cross2 1 , base 

b end 

ru 



T.r-V ■ — 

| J: j. pro draw_cross2_event , event 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
frame_no, file__name,$ 

total_f rame_no, time_resolution, x_size, y__size, box_size 
common markov2 , yes_no_values3 / added by pvc3 

stateholder5=WIDGET_INF0 (event .handler , /Child) 

WIDGET_CONTROL, stateholderS , GET_UVALUE=state 

WIDGET_CONTROL , state . draw, GET_VALUE=win_id2 

WIDGET_CONTROL , state . numberl , GET_VALUE =numbe r_l 

WIDGET_CONTR0L , state . number2 , GET_VALUE=number_2 

WIDGET_CONTROL, state . binsize , GET_VALUE=bin_size 

WIDGET_CONTROL , state . xmax, GET__VALUE=xmax 

WIDGET_C0NTROL , state . xmin, GET_VALUE=xmin 

IF ((number_l gt cell_no) OR (number_l It 1)) OR ( (number_2 gt cell_no) OR 
(number_2 It 1) ) THEN BEGIN 

mess=WIDGETJVIESSAGE (string (' Invalid cell number entered. You only 

chose ' , byte (cell_no) , 1 cells ' ) , /ERROR) 
ENDIF ELSE BEGIN 
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wset, win__id2 /changed by pvc3 

draw_cgram2 , yes_no_values3 , numberJL, number_2 , bin_size, cell_no, 
time_resolution, $ 

xmin, xmax 

ENDELSE 

end 

; Name: draw_raster 

; Description: This will draw spikes for those points that pass the test. 
; Working: It sees if the pixel vs time graph has any peaks and then simply 
draws a line 

; from the x-axis to represent an 'on' or spike. 

; It draws the plots one on top of another in a raster form. 

pro draw_raster 

common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 

f r ame_no , f i 1 e_name , $ 

total_frame_no, time_resolution, x_size, y_size, box_size 
common choose_cell_convert_delta, filepathl 

Q a=bytarr (f rame__no+l) 

H loadct f 27 

f% !p.multi=[0, 1, cell_no] 

$ window, 1, xsize=600, ysize=600 , title=strcompress (strmid (systime ( ) , 3, 8)) 

; Hack to print out raster values to file 
pi f ilename=pickf ile (/write, f ile= ' Location_Spikes_#_ ' ) 

CO if (filename eq ' 1 ) then begin 

55 MESS=WIDGET_MESSAGE ( 'Data will not be saved! ■ , / INFORMATION) 

}*&< cd, ' ' 

f ilename= ' rasterdata . txt ' 

f% endif 



■KM- 



[Z openw, outunit, filename, /GET_LUN 

H printf, outunit, FORMAT = '("Number of Frames: ",10," Number of Cells: 

*™ " , 10) 1 , f rame_no, cell_no 

printf, outunit, 1 1 

printf, outunit, 'Cell No No Spikes Spike Times (frame)' 

for cell_number = 0, cell_no - 1 do begin 

if (cell_number le 8) then begin ;; 8 because we add one to the 
cell number when we print 

printf, outunit, FORMAT = '("Cell ",I0,":",I," ",$)', 
cell_number+l, total (yes_no_values (cell_number , *) ) 
endif else begin 

printf, outunit, FORMAT = '("Cell ",I0,":",I," ",$)', 
cell_number+l, total (yes_no_values (cell_number , *) ) 
endelse 

; ; code to write out only those frames that have spikes 
if (total (yes_no_values (cell_number , *) ) ne 0) then begin 
;this only draws cells that have activity 

plot, a, xstyle=4, ystyle=4, ymargin= [0 , 0] , xmargin= [15 , 0] 
;xstyle and ystyle avoid display of x and y axis respec. see 2-8R 

xyouts, - (f rame_no/10 . 000) , 0, cell_number+l , charsize=0.8 
/this prints the cell numbers in front of the cells 
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for frames=0, frame_no-l do begin 

if (yes_no_values (cell_number, frames) ge 1) then begin 

; to read yes_no 

plots, [frames, frames], [.75, 0], color=3 
printf, outunit, FORMAT = '($,1)', frames 

endif 

endfor 

endif 

;; code to write yes_no_values instead- i.e. writes the whole binary 
spike train for each cell 

; printf , outunit, 1 1 

; printf , outunit, transpose (yes_no_values [cell_number , *] ) 
printf, outunit, 11 ; newline 

endfor 

CLOSE, outunit 
FREE_LUN, outunit 
!p .multi=0 

M; end — 

Q; NAME: draw_signif icance_raster . pro 

. 

\ DESCRIPTION: This program looks at the many_one and two_many significance 
%4 ^_ ^_ ^_ see 

ff ; whether many cells fire in one frame or a pair of cells fire together in 

'? many frames . 

■iS ; it then colors the two cases differently on the raster plot. 

10 ; This says nothing about the significance of the correlation. For that, 

% you must look at the 

U ; distribution curve or the data file that pops up after calling many_one or 

Id two_many 

P 
iff II 

IZ pro draw_signif icance_raster, yes_no_signif icance 

* A common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common with_create_dist , cells_active ; cells_active is the 

number of active cells in the movie 

DRAW THE BASE LINES FOR THE CELLS AND PLOT THE 



SPIKES ON THEM . 

a = bytarr(frame_no + 1) /plotting 'a' will make the line on 

which the spikes will be for each cell 
loadct, 27 

!p.multi=[0, 1, cells_active] ,-this plots cells_active many single lines 
one above another 

window, /free, xsize=600, ysize=600 , title= ' Significance Raster Plot' 
;the title is the date 

for cell_number = 0, cell_no - 1 do begin 

if (total (yes_no_signif icance (cell_number , *) ) ne 0) then begin 
;this only prints cells with activity 

plot, a, xstyle=4, ystyle=4, ymargin= [0 , 0] , xmargin= [15 , 0] 
;xstyle and ystyle avoid display of x and y axis respec. see 2-8R 
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xyouts, -(frame_no / 10.000), 0, cell_number + 1, color=2 
;this prints the cell numbers in front of the cells 
for frames = 0, framejno - 1 do begin 

COLOR THE SPIKES RED IF THEY HAPPEN TO FULFILL 

i _ , 

THE CRITERIA 

if (yes_no_signif icance (cell_number , frames) ge 1) then 

begin ;to read yes_no 

plots, [frames, frames], [.75, 0], 
color= (yes_no_signif icance (cell_number, frames) -l)*200+2 ;this makes the 

correlated spikes red and others grey 

endif 

endfor 

endif 
endfor 

PRINTING A DOTTED BLUE LINE OVER THE MANY CELLS THAT FIRE 



TOGETHER IN ONE FRAME_ 



IJ endfor 



for frames = 0, framejao - 1 do begin 

if total (yes_no_signif icance (*, frames) gt 1) gt 0 then begin 
plots, [frames, frames], [200,0], color=3 , linestyle=l 
endif 



!p .multi=0 



Name : draw_spikes 

Description: This will draw spikes for those points that pass the test. 
Working: It sees if the pixel vs time graph has any peaks and then simply 
*'* draws a line 

? ; from the x-axis to represent an 'on 1 or spike. 

H 

pro draw spikes, number_JL 

o 

RJ common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
P f r ame_no , f i 1 e_name , $ 

total_frame_no, time_resolution, x_size, y_size, box_size 

common share_w_j>ixel_vs_time_widget_com, y_min, y_max, spikejnin, spike_max, 
x__max, x_min 

!p .multi=0 

plot, f indgen(frame_no) *time_resolution, pixel_array (number_l-l, *) , 

yrange= [y_min, y__max] , $ 

xrange= [x_min, x_max] , title=number_l , ytitle='% Change in Index', 

xtitle='Time (days)', FONT=-l 

for frames=l, frame_no-l do begin ;this doesn't consider a spike 

at the 1st frame 

if (yes_no_values (number_l-l, frames) ne 0) then begin 

plots, [frames*time_resolution, f rames*time_resolution] , 
[spike_max, spike_min] 
endif 

endfor 

end 



; Name: draw_spikes_many_cell_delta 
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; Synopsis : 

; Description: This program plots the deltaF/F values of a single cell per time. 

pro draw_spikes_many_cell_delta, delta_f __valu.es , cell_number , 
f rame_number , time_resolution 

!p.multi=0 

window, 2, title= ' Superimposed Plots' 

plot, findgen(frame_number) *time_resolution / delta_f_values (0 , *) , 
yrange= [min (delta_f_values) , max (del ta_f_va lues) ] , $ 

title='All Stocks', ytitle= ' % Change in Index', xtitle='Time (days)' 



xyouts, -6, 0, 1 

if (cell_number gt 1) then begin 

for cell_numb=l, cell_number-l do begin 

oplot, findgen(frame_number) *time_resolution, 
(delta_f_yalues (cell_numb, *) ) 

a= (delta_f_values (cell_numb # *) eq 
max( (delta_f_values (cell_numb, *) ) ) ) 
a=where(a eq 1) 

xyouts, a,max( (delta_f_values (cell_numb, *) ) ) , cell_numb+l 

r ;" : endfor 

™- endif 

m 

%v# end . 



W ; Name: draw__spikes_many__pages 

III 

\|pro draw_spikesjnany__pages, pixel_array, yes_no_values , cell_no, 
f'fli f rame_number, time_resolution 

fh 

^ common share__w_pixel_vs_time_widget_com, y_min, y_max, spike_min, spike_max 

jj sum_array=intarr(cell_no) ;find which cells are 

M active 

nj for counter=0, cell_no-l do begin 

p sum_array (counter) =total (yes_no__values (counter, *) ) 

M> endfor 

active_cells=total(sum_array gt 0) ;this has the location of the 

cells that are active 

if active_cells eq 0 then begin 

mess=WIDGET_MESSAGE ( ' There are no active stocks', /error) 
endif else begin 

active_cell_subscripts=where (sum_array gt 0) ;this is equal 

to the number of active cells 

number = (active_cells / 9) + 1 ; number 

of windows that have to be opened 

for t = 0, number - 1 do begin ;open the 

windows 

window, t, xsize=500, ysize=500 

endfor 

!p.multi= [0,3,3] ;it' 11 print three rows of 

three cells 

for index=0, active_cells-l do begin 
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wset, (index) / 9 

cell_number=active_cell_subscripts (index) 
plot, f indgen(frame_number) *time_resolution, 
pixel array (cell_number , *) , yrange= [y_min, yjnax] , $ 

title=cell_number+l, ytitle='% Change in Index', xtitle='Time 

(days) ' 

for frames=0, f rame_number-l do begin 

if (yes_no_yalues (cell_number, frames) ge 1) then begin 

plots, [frames*time_resolution, f rames*time_resolution] , 

[ spike_max , spike_min] 
endif 

endf or 

endf or 
endelse 

"p.multi=0 ; this returns plotting to 

systems default 

end m _ 

; Name: f ilter_median_subtractive 

u ;; Parameters: takes an input array (the array of cells x frames dF/F values), a 
window size for 

y ; filtering, and an output array (an empty variable which will be the 

W filtered input passed back) . 

Its:?-- ; 

IU ; 

Si 

flj pro filter_median__subtractive, input_array, median_window_size, output_array 

IB 

output_array = f loat (input__array) 
!\ num_cells = n_elements (input_array (* , 0) ) 

num_f rames = n_elements (input_array (0 , *) ) 

temp array = f ltarr (num_f rames + (2 * median_window_size) ) 
IU for i = 0, (num_cells - 1) do begin 

Q temp_array (median_window_size:median_window_size + num_f rames - 1) = 

hi reform ( input_array ( i , * ) ) , 

temp__array = (temp_array - median (temp_ar ray , median__window_size) ) 
output__array ( i , * ) = temp_array (median_window_size :median_window_size 
+ num_f rames - 1) 
endf or 

end . — 

/ Name: f ind_distance 

; Description: This procedure is used by draw_cgram to make cross correlograms 
and auto correlograms. 

; Timearrayl and timearray2 are arrays which hold the times where the spikes 
occur. The output 

; of this procedure is an array of the difference of the second spike tram from 
the first cell's. 

; It does this by taking the spike of the first cell eg 13 54 66 154 and looks 
at the second spike 

train's timearray which may look like: 14 55 67 155 
outputarray (0, 0) =14-13=1 
outputarray (0, 1) =14-54=-40 
outputarray (0, 2) =14-66=-52 
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outputarray (0, 3 ) =14 -154= -14 0 
outputarray (1, 0) =55-13=42 
output array (1, 1) =55-54 =1 
and so on 

this will show a peak at 1 in the cross correlogram for 1 to 2 . 
function f ind_distance, timearrayl, timearray2 

;this function takes two files and finds how their elements differ by 
subtracting all of the elements 

; of one from the other. Output is a file with a measure of their deviation from 
each other 

sizel=n_elements (timearrayl) 
size2=n_elements (timearray2) 
outputarray=intarr (sizel, size2) 
index=0 

for i=0, sizel-1 do begin ;one by one subtracting the elements 

of filel 

;from the whole of file2 
output array (index, *) =timearray2 -timearrayl (i) 
I , index=temporary (index) +1 

*~ endfor 

P ;for j=0, sizel-1 do begin 

;for t=0, size2-l do begin 
U ;outputarray ( j , t) =timarray2 (t) -timearrayl ( j ) 

ft! ; endfor 

\| ; endfor 

|;0 return, outputarray 
m end 

:ii 

U 

u 

•».;*• 

ru 



NAME: f ind_matches .pro 
INTRODUCTION : 

This calculates the correlation coefficient of two cells with a binsize 



INPUTS : 

celll and cell2 are the unidimensional arrays of the timearrays of the 
cells . 

winsize is the size of the window. THIS MUST BE AN ODD NUMBER. 

DESCRIPTION; 

The c coeff is calculated as follows. 

1. The binary spikes are converted into time arrays where the values 

of the 

array correspond to the indices where the spikes occurred. 

For example, if celll=[0, 0, 0, 0, l, 0, 0, 1, 0, 0, l, l] 
timel=where (celll gt 0) 
IDL> print, where (celll gt 0) 

4 7 10 11 

2. The IDL library function strpos is used to find whether any of 
the time array 

; values are the same. If there are two time array values that are the same, 
that means that both of 

; the cells we are analyzing had spikes at the same time. That counts as a match 
in the formula to 

; calculate correlation coefficients: 

c_coeff= number of matches (A- >B) / total spikes of A 
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; 3. f indjnatches counts the number of times there are congruent spike 

positions. The 

number of these occurences = matches 

4. The total number of spikes in cell 1 =lengthl 

5 Cross correlation coef =matches/lengthl 

function f ind_matches , celll, cell2, winsize 

matches=0 ; matches is the number of coincident spikes 

timel=where (celll gt 0) ;gives the location of where celll is 

greater than 0 

time2=where (cell2 gt 0) ;gt means greater than. similarly->ge is 

greater than equal 

lengthl=total (timel gt 0) /lengthl gives the number of elements in the 



array timel 
1^ winsize) as a spike 



/which is also the number of spikes in celll 
;the same window (whose size is defined by 

;from cell2. It is calculated in the module below 



; FIND THE NUMBER OF COINCIDENT 

IM SPIKES 

M 

CO for 1=0, lengthl -1 do begin 

HI window_size= (winsize-1) /2 ;a window is created x before a spike and x 

s after. Therefore making it 2x+l large 
| ?! : ; lower_window_limit=timel (i) -window_size 

upper_window_limit=timel (i) +window_size 
ft% count=0 /count is the number of times a spike from 

Jl"? celll occurs in 

! .j for m=lower_window_limit , upper_window_limit do begin 

y j temp=strpos (time2 , m) 

M ; templ= (temp) +1 /there is a zero if there is a match 

therefore we have to add 1 

if (total (tempi) ne 0) then count=count+l /in this case there is a 
coincident spike 
endf or 

if (count gt 0) then matches=matches+l 

endfor 

; CALCULATE THE CORRELATION 

COEFFICIENT 



coef f=matches/ (lengthl) /lengthl is always greater than 0 as we only 

do this for active cells 
return, [coeff] 
end 

; gen_sig_widget 
pro gen_sig_widget 

gen_sig_base=WTDGET_BASE (/ COLUMN, title= 1 General Significance' ) 
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gen_sig_num_trials=CW_FIELD(gen_sig_base, /RETURN_E VENTS , / INTEGER , 
TITLE= 1 Number of random trials to run: 1 , VALUE=1000, UVALUE=0) 

gen_sig_window_size=CW_FIELD(gen_sigjDase, / RE TURN_E VENTS , / INTEGER , 

TITLE= ' Window size:', VALUE=0, UVALUE=0) 

gen_sig_button=WIDGET_BUTTON (gen_sig_base , VALUE= » Find significance 1 

UVALUE=2) 



WIDGET_ CONTROL , /realize, gen_sig_base 

gen s ig_s tat e={gen_sig_num_t rials :gen_sig_num_t rials , $ ^ 

gen_sig_window_size :gen_sig_window_size, $ 

gen_s ig_but ton : gen_s ig_but ton } 
WIDGET_CONTROL, WIDGET_INFO (gen_sig_Jbase , /CHILD) , 

SET_UVALUE=gen_sig_state 

xmanager, ' gen_sig_widget 1 , gen_sig_base 

end ____ 

; Name: gen_sig_widget_event 



pro gen_sig_widget_event , event 

u : common mother_com, pixel_array, yes_no_values, coef, location, cell_no, 

f rame_no , f i le_name , $ 
2 total_frame_no, time_resolution, xjsize, y_size, box_size 

p common with_create_dist , cells_active 

W gen_sig_stateholder=WIDGET_INFO (event . handler , / CHILD) 

;?; J WIDGET_CONTROL, gen_sig__stateholder , GET_UVALUE=gen_sig_state 

II WIDGET_CONTROL , gen_sig_state . gen_sig_num_t rials , GE T_VALUE = no_o f _t i me s 

fS WIDGET CONTROL , gen sig state . gen_sig_window_jsize , GET_VALUE=window_size 

i\f,i< — — — 



if (event, id eq gen_sig_state .gen_sig_button) then begin 

if ( (no_of_times le 1) or (window_size le -1)) then begin 

error_dialog = WIDGET_MES SAGE (» Invalid fields specified! ' , 

;; /error) 
HI endif else begin 

Q ;; We write all of our data into the following file in the 

|»*< function f ind_general_p 

filename3 = " GeneralJ3tatistical_Data 1 

filename3 = pickfile (/write, f ile= 1 General_Stats_# ' ) 

if (filename3 eq ■ 1 ) then begin 

mess = WIDGET__MESSAGE ( 1 The Data will not be saved.', 

/ INFORMATION) 

filename3 = 1 General_Statistical_Data ■ 

endif 



WIDGET_CONTROL, /HOURGLASS 

; ; creating the random distribution 

times_repeated = 2 

no_spikes = total (yes_no_values) 

distribution_array = intarr (no_of_times) ; ; array of hits for 

each random trial 

for t = 0, no_of_times - 1 do begin 

make_random_data, seed, random_array, num_active_cells 
distribution_array (t) = 
count_random_matches (random_array , times_repeated, window_size) 
endf or 
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; ; count the number of matches in the real data now 
true_matches = 0 

if (window_size eq 0) then begin 

for frame_counter = 0, frame_no - 1 do begin 

matches = total (yes_no_values (* , f rame_counter) ) 
;this gives us the n in our combination formula above 

if (matches ge times_repeated) then begin 
matches_pair = factorial (matches) / 
(factorial (times_repeated) * factorial (matches - times_repeated) ) ; finds 
number of possible pairs 

endif else begin 

matches_pair = 0 ; cant find 

factorial of a negative number 

endelse 

true_matches = true__matches + matchesjjair 
; increases the matches count for the true data by that of this frame 

endfor 
endif else begin 

;; first compress spike trains according to window size 
binned_yn_vals_array = intarr (cell_no, frame_no - 



(window size * 2)) 
y +1) do begin 



for frame_counter = window_size, frame_no - (window_size 

for window_counter = -window_size , window_size do 
4 :f; | begin 

ks for cell_counter = 0, cell__no - 1 do begin 

;1 if (yes_no_yalues [cell_counter, 

m frame counter + window_counter] ge 1) then begin 

%M — 

m binned_yn_vals_array [cell_counter , f rame_counter - window_size] = 1 

f endif 
I"* endfor 
H : endfor 
O endfor 

fjj ; ; then look for matches in the binned spike trains 

p (although a single spike can fall into many bins) 
yt. for f rame_counter = 0, 

n elements (binned_yn_vals_array [0, *] ) - 1 do begin 

matches = total (binned_yn_vals_array (* , 
f rame_counter) ) ; finds number of hits at this frame 

if (matches ge times_repeated) then begin 
matchesjpair = factorial (matches) / 
(factorial (times_repeated) * factorial (matches - times_repeated) ) ; finds 
number of possible pairs 

endif else begin 

matchesjpair = 0 ; cant find 

factorial of a negative number 

endelse 

true_matches - true_matches + matches J>air 
/increases the matches count for the true data by that of this frame 

endfor 
endelse 



; ; determine connectivity stats 
distances_array = [0.0] 
total connections = 0 
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connect ions_ar ray = intarr (cell_no) 
if (window_size eq 0) then begin 

for cell_l = 0, cell_no - 2 do begin 

for cell_2 = cell_l + 1, cell_no - 1 do begin 
temp_cell = yes_no_values (cell_l , *) * 

yes no values (cell_2 , *) 

~ ~ if ( (total (temp_cell) ) ge 1) then begin 

connections_array [cell_l] = 

connect ions__array [eel 1_1] + 1 

connections_array [cell__2] = 

connect ions_array [cell_2] + 1 

total_connections = total_connections 
+ 1 ; ; do this so a connection is only counted once, if it matters 

cojmection_distan.ee = 
sqrt (((double (location (cell_l) .coord [0] ) - double (location (cell_2) . coord [0] ) ) * 
2) + ( (double (location (cellJL) .coord [1] ) - double (location (eel 1^2) . coord [1] ) ) 
2)) 

distances_array = 

[distances_array, connection_distance] 

endif 

endf or 

s . endf or 

endif else begin 
U for cell_l = 0, cell_no - 2 do begin 

P for cell_2 = cell_l + 1, cell_no - 1 do begin 

P temp_cell = binned_yn_vals_array (cell_l , *) 

0 J * binned yn vals_array (cell_2 , *) 

y ~ " if ( (total (temp_cell) ) ge 1) then begin 

ffi, connect ions_array [eel 1__1] = 

m connections_array [cell_l] + 1 

1:V " connections_array [cell_2] = 

$ 

connections_array [cell_2] + 1 
I"* total_connections = total_connections 

t"! + 1 ;; do this so a connection is only counted once, if it matters 

P ' connection_distance = 

fij sqrt {((double (location (cell_l) .coord [0] ) - double (location (cell_2 ). coord [0] ) ) * 

P 2) + ( (double (location (cell_l) .coord [1] ) - double (location (cell_2 ). coord [1] ) ) 

U 2 ) ) 

distances__array = 

[distances_array / connect ion_di stance] 

endif 

endf or 

endf or 
endelse 

/ ; Begin printing the data 
close, 1 

openw, 1, filename 3 

printf, 1, 1 STATISTICAL DATA ' 

printf, 1, ' ' 

printf, 1, 'Random Test Information' 

printf, 1, ' ' 

printf, 1, 'How many times a pair is together: ' , 
strcompress (times_repeated) 

printf, 1, 'Number of iterations in tests: 
strcompress (no_of_times) 
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printf, 1, 'Window size: 
strcompress (window_size) 

printf, 1, 'Total number of active stocks: " , 

strcompress (round (cell s_active) ) 
printf, 1, ' 1 

if (max(distribution_array) eq 0) then begin 

mess = WIDGET_MESSAGE ( ' Two cells never fire together 

randomly! Try the MANY CELLS AT ONCE option looking at 2 cells firing together 

with more than 100 iterations.', /error) 

printf, 1, 'ERROR!' 

printf, 1, 'Two cells never fire together randomly. 
Nothing can be said about the significance of this data' 

real J = -1 
endif else begin 
loadct, 2 7 

hist = histogram (distribution_array, binsize=l, min=0, 
max= (max( [truejnatches , distribution_array] ) + 2) ) 

window, /free, title= 'Distribution for General 

Significance' , 

plot, hist, xtitle=' number of hits', yt it le= • frequency 

plots, [true_matches, true_matches] , [0, 

max ( [1, total (distribution_array eq true_matches) ] ) ] , color=12 ;this draws a 

W line where the actual number of matches lies 

C3 

p spikes_per_cell_j>er_second = (no_spikes / cells_active) 

flj / (frame_no * time_resolution) 

m 



gen sdev) 

ifi — 



stats = moment (distribution_array, sdev=gen_sdev) 
p_value = gaussjpdf ( (true_matches - stats (0)) / 

realjp = 1 - p_value 

printf, 1, 'General Statistical Information 

printf, 1, 1 

printf, 1, 'Total number of frames: 



Qi strcompress (frame_no) 



printf, 1, 'Time between frames: 
strcompress (time_resolution) 

printf, 1, 'Total number of stocks: 



strcompress (cell_no) 
strcompress (no_spikes) 
strcompress (stats (0) ) 
strcompress (stats (1) ) 
strcompress (gen_sdev) 



printf, 1, 'Total number of spikes: 

printf, 1, 'Mean expected matches: 

printf, 1, 'Variance: 

printf, 1, 'Standard deviation: 



printf, 1, 'Actual matches: 
strcompress ( true_matches) 

printf, 1, 'Actual/expected: 
strcompress (true_matches / stats (0)) 

printf, 1, 'Standard error for ratio: 

strcompress (gen_sdev / stats (0)) 

printf, 1, 'Mean spikes per stock: 
strcompress (no_spikes / cells_active) 
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printf, 1/ 'Spike firing rate: 
strcompress (spikes_per_cell_per_second) 

printf, 1, 'Significance p-value: 

strcompress (real _j>) 

printf, 1, ' ( IMPORTANT ! If this is exactly .5, it could 

be actually very significant so look at actual/expected!) ' 

printf, 1, ' ' 



printf, 1, 'Stock Connectivity Information' 

printf, 1, ' ' 

printf, 1, 'Total Number of Connections: 

strcompress (total_connections) 

connect ed_cells = where (connect ions_array) 
num_connected_cells = n_elements (connected_cells) 
if (n_elements (connected_cells) ge 2) then begin 
stats_conn_only = 
moment ( connect ions_array (connected_cells) , sdev=sdev_conn_only) 

mean_connections = stats_conn_only [0] 
normalized_slice_connectivity = (mean_connections 

/ num connected_cells) 

printf, 1, 'Mean Number of Connections per Stock: 
j ' , strcompress (mean_connections) 

P printf, 1, 'Standard Deviation: 

P ' , strcompress (sdev_conn_only) 

flj endif else begin 

\| normalized_slice_connectivity = 0 . 0 ; ; no 

fn connected cells! Can't have just one connected cell. 
fi% endelse 

printf, 1, 'Normalized Connectivity: 
s t r compr e s s ( norma 1 i z ed_s lie e__conne c t i vi ty ) 

printf, 1, 'Number of Silent Stocks: 
strcompress (cell_no - num_connected_cells) 

printf, 1, 'Number of Coactive Stocks: 
strcompress (njslements (connected_cells) ) 

printf, 1, 'Minimum Number of Connections: 
strcompress (min (connect ions_array) ) 

printf, 1, 'Maximum Number of Connections: 
strcompress (max (connect ions_array) ) 

if (n_elements (distances_array) gt 2) then begin 
distance_stats = 
moment (distances_array [1 :n_elements (distances_array) - 1] , 
sdev=sdev_distances_array) 

printf, 1, 'Mean Connection Distance: 
1 , strcompress (distance_stats [0] ) 

printf, 1, 'Standard Deviation: 
? , strcompress (sdev_di stance s_ar ray) 
endif 
endelse 

; ; Close the file we are writing, reopen it if small enough 
close, 1 

openr, 1, filename3 
stat_stuff = fstat(l) 
file_size = stat_stuf f . size 
close, 1 



m 
o 
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if (file_size gt 8112) then begin 

mess=WIDGETJV!ESSAGE ( 'File is too large to display 
through a widget. Open it manually. If you did not save it, it is named 
General_Statistical_Data' , /INFORMATION) 

endif else begin 

xdisplayfile, filenames, title = "General Statistical 
Data", group = event. top, width = 75, height =50 
endelse 
endelse 

endif 

end 



This file reads MS Excel file with stock data of the following format 

date tickerl ticker2 

value closel close2 

This file creates the following arrays: 

symbol_array (string) - first row of the excle file minus first value 
date_array (long or string or date) - first column of the excel file minus 
the first value 

;; pixel_array (float) - all the rest 

; ; pixel_array later gets transformed by calculating deltaF/F - under the same 
M: name 

D 

■m pro load and convert excelfile 

TflS' 

HI 

\l common mother_com, pixel_array, yes_no_values, coef , location, cell no, 

* frame_no, file_name,$ ~ 

total_frame_no, time_resolution, x_size, y_size, box_size 
common test, str_ing, state3 
common old_skool_data, original_data 
\A common flags, cells_def ined, spikes_def ined, single_plot defined, 
j,* correl__coef_de fined 

O common with_choose_cells_com, thresholdl, threshold2 , rms_threshold 
fU common textfile_vars, text_flag, diode_array, max_num_of_diodes 
J,;;j. common stockdata, symbol_array, date_array 

message_dialog=WIDGET_MESSAGE(<'This procedure loads stock data, it creates 
symbol_array, date_array, and pixel_array " , /INFORMATION) 

filename=dialog_jpickfile(/read, f ile= (' stocks . slk ») , getj>ath=f ilepathl) 
if (filename eq ") then begin 

message_dialog=WIDGET_MESSAGE("No data read.", / INFORMATION) 
endif else begin 
close, /all 
openr, 1, filename 



;use excel import function to fill all the required arrays 
symbol_array= read_sylk (filename , / ARRAY , nrows = 1 , startcol=l) 
date_array=read_sylk (filename, /ARRAY, ncols=l, 
startrow=l, /uselongs) 

data = read_sylk( filename, / ARRAY, startrow=l, startcol=l) 
close, /all 

cell_no = n_elements (symbol_array) 



NY02:354955.1 



-55- 



f rame_no - n__elements (date_array) 

WIDGET_CONTROL , state3.name, SET_VALUE=strmid (filename, 
strlen (f ilepathl) ) 

WIDGET_CONTROL, state3 . frame, S E T_VALUE = f r ame_no 
WIDGET_CONTROL , state3 . time, SET_VALUE=time__resolution 

;; build 1 pixel_array 1 

pixel_array = f ltarr (cell_no, frame_no) 
for j=0, frame_no -1 do begin 

for i = 0, cell_no - 1 do begin 

pixel_array [i, j ] = data[j,i] 

endfor 

endfor 

original_data = pixel_array 

;the following finds delta F over F 

deltaf =f ltarr (cell_no, f rame_no) 

for i = 0, cell_no - 1 do begin 

for j=l, frame_no-l do begin 

deltaf [i, j] =100* (pixel_array [i , j ] - pixel_array [i , j - 

1] ) /pixel_array [i, j -1] 
m endfor 
g endfor 

p for i = 0, cell_no - 1 do begin 

deltaf [i, 0] =0 

B endfor 

m 
m 

* ; ; initialize other variables 

H yes_no_values = intarr (cell_no, frame_no) 

!#: coef = f ltarr (cell_no, cell_no) 

p rms_threshold = f ltarr (cell_no) 

I'jij thresholdl = f ltarr (cell_no) 

threshold2 = f ltarr (cell_no) 
for cell = 0, cell_no - 1 do begin 
rms_threshold [cell] =2.0 
thresholdl [cell] =2.0 
threshold2 [cell] =3.0 

endfor 

cells_def ined = 1 
spikes_def ined = 0 
single_plot_def ined = 0 
correl_coef_def ined = 0 
box_size = 2 
x_size = 350 
y_size = 350 

; ; locations are assigned deterministically around a circle 
location = replicate ( {struct , coord : intarr (2) , size:0, 
half_side : 0 . 00} , cell_no) 

for i = 0, (cell_no - 1) do begin 
location (i) .size = cell_no 
location (i) .half_side = 1 
location (i) .coord [0] = cell_no * 5 + 1 
location (i) . coord [1] = cell no * 5 + 1 



pixel_array=deltaf 
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endf or 

text_flag = 1 

message_dialog=WIDGET_MESSAGE ( "Finished reading data . " , 
/ INFORMATION) 
endelse 

end 

; Name : make_binary 

function make_binary, delta_values , frame_no, cell_no, thresholdl, threshold2 

binary_values = intarr (cell_no, frame_no) 

for cell = 0, cell_no - 1 do begin 
Case 1 of 

(thresholdl [cell] gt 0): begin ; f or increasing 

spikes, the threshold will be greater than 0 

for frame=l, frame_no-l do begin ;it doesn't consider a 
spike at the first frame 

if (frame It frame_no-2) then begin ;for all of 

the peaks upto the last frame 

if ( ( (delta_values (cell, frame+1) ) - 
U (delta_values (cell, frame)) ge thresholdl [cell] ) OR $ 

13 (delta__values (cell, frame+2) -delta_values (cell, 

Q frame) ge threshold2 [cell] ) ) then begin 

f!J binary_values (cell, f rame) =1 

%3 end if 

fp endif else begin ;this takes care of a peak at last 

frame 

^ if ( (delta_values (cell, f rame_no~l) ) - 

f (delta_values (cell , frame_no-2) ) ge thresholdl [cell] ) then begin 

I"* binary_values (cell, f rame_no-2) =1 

M endif 

0 endelse 

jffj endfor 

fli end 

u 

(thresholdl [cell] It 0): begin ;for decreasing 

spikes, the threshold will be less than 0 

for frame=l, frame_no-l do begin ; it doesn't consider a 
spike at the first frame 

if (frame It frame_no-2) then begin 

if ( ( (delta_values (cell, frame+1) ) - 
(delta_values (cell, frame)) le thresholdl [cell] ) OR $ 

(delta_values (cell, frame+2) -delta_values (cell, 
frame) le threshold2 [cell] ) ) then begin 

binary_values (cell, frame) =1 

endif 

endif else begin ;this takes care of a peak at last 

frame 

if ( (delta_values (cell, f rame_no-l) ) - 
(delta_values (cell, f rame_no-2) ) le thresholdl [cell] ) then begin 

binary_values (cell, frame_no-2) =1 

endif 
endelse 

endfor 
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end 

endcase 

endfor 



return , binary_yalues 
end 

; Name : make_random_data . pro 

; Description: This program looks at the input cell's activity and creates a 
random set of data 

; based on this data, by rotating each cell's spike train by a random amount. 
; Please note that this random array contains ONLY those cells which spike 

at least once! 



pro make_random_data, my_seed, random_array, num_active_cells 

common mother_com, pixel_array, yes_no_values, coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_re solution, x_size, y_size, box_size 
common with_create_dist , cells_active 

/find the number of spikes each cell of the true data has 
no_spikes = intarr (cell_no) ;this is the 

array with the number of spikes for each cell 



W for cell = 0, cell_no - 1 do begin 

C3< no_spikes (cell) = total (yes_no_values (cell , *) ) 

B endfor 

m 

Si cells_active = total (no_spikes gt 0) 

|pj num_active_cells = cells_active 

\k random_array = intarr (cell s_active, frame_no) 

random cell index = 0 

for cell_counter = 0, cell_no - 1 do begin 

if (no_spikes (cell_counter) ge 1) then begin 
|3 ; ; pick a random int between 0 and f rame_no 

fjj random_time_shif t = f ix (randomu (my_seed) * frame_no) 

p random_array [random_cell_index, *] = 



shift (yes_no_yalues [cell_counter, *] , random_time_shif t) 

random__cell_index = random_cell_index + 1 

endif 

endfor 

end 



; Name : 

; make_single_binary 
; Description: 

; This program comes under the single plot call. It is called when the user 

enters 2 thresholds that 

; are different from the 2 thresholds that calculated the spikes for the whole 
movie (under the choose threshold call) . 

; Make_binary calculated the spikes for the whole movie while this function is 
called make_single_binary . 

; This calculates the spikes only when told to do so for a single cell under the 
single plot widget. Therefore it is called 
; make_single_binary . 
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It goes through the deltaF/F values (delta__values) and sees if the 
absolute increase in delta_value between adjacent frames 

/ is greater than or equal to the threshold for positive spikes and less than or 

equal to the threshold for negative spikes . 

; We do not consider spikes at the first frame. 

/ 

; Explanation of variables: 

/ delta_values : the array with the actual delta F/F values 

; yes_no_values : the converted binary array. 

cell_number: the cell whose spikes are being found. 
; frame_no: total number of frames 

/ threshl: threshold between two adjacent frames 

; thresh2 : threshold between three adjacent frames 



function make_single_binary , delta_values, yes_no_values , cell_number, frame_no, 
threshl, thresh2 

yes_no_values (cell_number-l , *) =0 
Case 1 of 

(threshl gt 0) : begin ;for increasing spikes, the threshold 

H : will be greater than 0 

O for frame=l, frame_no-l do begin ;this doesn't consider a spike at the 

Q first frame 

p if (frame It frame_no-2) then begin 

if ( ( (delta_yalues (cell_number-l , frame+1) ) - 
(delta_values (cell_number-l , frame)) ge threshl) OR $ 
^~ (delta_values (cell__number-l , frame+2) - 

fl delta_values (cell_number-l, frame) ge thresh2)) then begin 
W yes_no_values (cell_number-l, frame) =1 

a endif 

H ; endif else begin ;this takes care of a peak at point 1 

if ( (delta_values (cell_number-l , frame_no-l) ) - 
p (delta_values (cell_number-l, frame_no-2) ) ge threshl) then begin 
llj yes_no_values (cell_number-l, frame_no-2) =1 

f?s endif 
! ; ? endelse 
endfor 

end 



(threshl It 0) : begin ; f or decreasing spikes, the threshold 

will be less than 0 

for frame=l, frame_no-l do begin ;this doesn't consider a spike at the 
first frame 

if (frame It frame_no-2) then begin 

if ( ( (delta_values (cell_number-l , frame+1) ) - 
(delta_values (cell_number-l, frame)) le threshl) OR $ 

(delta_yalues (cell_number-l, frame+2) - 
delta_values (cell_number-l, frame) le thresh2) ) then begin 

yes_no_values (cell_number-l , frame) =1 

endif 

endif else begin ;this takes care of a peak at point 1 
if ( (delta_values (cell_number-l, f rame_no~l) ) - 
(delta_values (cell_number- 1 , f rame_no-2) ) le threshl) then begin 

yes_no_values (cell_number-l , frame_no-2) =1 

endif 
endelse 
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endf or 

end 

endcase 

return, yes_no__values 
end 



many_c ell s_one_f rame 

Creates a random distribution (via count_random_matches_many and 
make_r andom_da t a ) 

to compare the real data to- counts the number of times a minimum number of 
cells fire in one frame, or a set of frames given by the window size. 

function many_cells_one_f rame , synchronous_cells , no_iterations , ge_or_eq_test , 
window__size 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
frame_no, f ile_name, $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common with_create_dist , cells_active 
common withjnany, filename 

common stockdata, symbol_array , date_array 

/common with_correl_map_plan, cluster_array ;uncomment this if you 

& want networks to have different linestyles in the correl_map _plane . . 

P ; ; Initialize variables 

jlj connections_array = bytarr (cell_no, cell_no) ; if two cells are 

\j significantly connected, this is 1; if not, 0 

ffi ye s_no_signif icance = yes_no_values ; the new colored raster 

fi plot is made from ye s_no_signif icance . This sets it to the original data 
*- cluster_array = intarr (cell_no, cell_no) ; this holds the frame in 

jf_ which the cells fire together 

!«* true_hits = 0 ; number of networks 

H 1 in the real data 

0 distances_array = [0.0] ; list of all network 

fyj connection distances 
pi areas_array = [0.0] 

L.= ; network areas 

size_array = intarr (frame_no) 
those networks ' synchronous_cells ' + in size 

cluster_no = 0 
first network in solid lines 



list of all 
this counts only 
this makes the 



; ; Creating the random distribution 

no_spikes = total (ye s_no_values) /this give the 

total number of spikes by cells in the slice 

random_distribution = intarr (no_iterat ions) /this array 

holds one random number of matches for every iteration 
for t = 0, no_iterations - 1 do begin 

make_random_data, seed, random^array, num_active_cells 
random_distribution (t) = count_random_jnatches_many (random_array, 
sy nc hr onous_c ells, ge_p r__e q_t est, wi ndow_s i z e ) 
endf or 

/ / Printing some of the data 
close, 1 
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0 



filename = pickfile (/write, file= (strcompress (string (synchronous_cells) ) + 
'_Stocks_Once ' ) ) 

if (filename eq ,! ) then begin 

mess = WIDGET__MESSAGE ( 1 This data will not be saved r , /INFORMATION) 
f i 1 ename = 1 Many_One . da t ' 

endif 

openw, 1, filename ;D0 THIS-this is the temp file used for 

writing into the xdisplay. It can be saved to a further name, 
printf, 1, 'Statistical data for the Many/Once test.' 

printf, 1, 'One hit is when 1 , strcompress (synchronous_cells) , 1 stocks 
fire in one frame' 

printf, 1, 'Number of iterations: ', no_iterations 

printf, 1, 'Window size: ', window_size 

printf, 1, 'Below is the list of stocks which spike in a given frame:' 

; ; Finding the real number of matches- first bins spikes according to 
window size, no binning done when window is 0 

binned_yn_vals_array = intarr (cell_no, frame_no - (window__size * 2)) 
for f rame_counter = window_size, frame_no - (window_size + 1) do begin 
for window_counter = -window_size , window_size do begin 
for cell__counter = 0, cell_jio - 1 do begin 

if (yes_no_values [cell_counter , frame_counter + 
window_counter] ge 1) then begin 

binned_yn_vals_array [cell_counter , f rame_counter - 

p window_size] = 1 
P endif 
f]j endfor 
S.J endfor 
endfor 

* yi ; ; Now looks for matches in the binned data- there may be false positives 

f for windows 

for f rame_counter = 0, n_elements (binned_yn_vals_array [0 , *] ) - 1 do begin 
H ; matches = total (binned_yn_vals_array (* , f rame_counter) ) 

|5 if ( ( (ge_or_eq_test eq 0) and (matches ge synchronous_cells) ) or 

Ilj ( (ge_or_eq_test eq 1) and (matches eq synchronous_cells) ) ) then begin 
□ true_hits = true_hits + 1 

|^ ( size__array [f rame_counter] = matches 

printf, 1, 1 1 

if (window_size eq 0) then begin 

printf, 1, strcompress (round (matches) ) , ' Stocks Spike 
in Frame: ', strcompress (frame_counter) 

correlated_cell_array = where (ye s — no_values (* , 

frame_counter) eq 1) 

endif else begin 

printf, 1, strcompress (round (matches) ) , ' Stocks Spike 
at Approximately Frame: 1 , strcompress (frame_counter) 

correlated_cell_array = where (binned_yn_vals_array (* , 

frame_counter) eq 1) 

endelse 

network_xcors = intarr (n_elements (correlated_cell_array) ) 
network__ycors = intarr (n_elements (correlated__cell_array) ) 
for index = 0 , matches - 1 do begin 

printf, 1, 1 Stock number: ' , 

correlated_cell_array (index) + 1 

loc = correlated_cell_array (index) 
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for window_counter = f rame_counter , f rame_counter + (2 * 

window_size) do begin 

if (yes_no_values (correlated_cell__array (index) , 
window_counter) ge 1) then begin 

yes_no_signif icance (correlated_cell_array (index) , window_counter) = 
yes__no_signif icance (correlated_cell_array ( index) , window_counter) + 1 

endif 

endfor 

network^xcors [index] = location (loc) . coord [0] 
network_ycors [index] = location (loc) . coord [1] 
printf, 1, ' X and Y: ' , location (loc) . coord 

printf, 1, 1 1 

for index2 = index, matches - 1 do begin 
if (index ne index2) then begin 
connect ion_di stance = 
sqrt ( ( (double (location (correlated_cell_array [index] ) . coord [0] ) - 
double (location (correlated_cell_array [index2] ) .coord [0] ) ) 2) + 
( (double (location (correlated_cell__array [index] ) . coord [1] ) - 
double (location (correlated_cell_array [index2] ) .coord [1] ) ) A 2)) 
N ! distances_array = 

U [distances_array / connection_distance] 
Q " ~ endif 

Q loc2 = correlated_cell_array (index2) 

jfjj connect ions_arr ay (loc, loc2) = 1 

cluster_array (loc, loc2) = cluster_no 

endfor 

endfor 



m 

f ; ; find area of many cells once network 

N ! if ( (n_elements (uniq (network_xcors , sort (network_xcors) ) ) ge 

M 3) and (n_elements (uniq (network_ycors , sort (network_ycors) ) ) ge 3)) then begin 
Q triangulate, network_xcors, network_ycors , triangles, 

T\,$ convex_hull ; ; here we find the convex hull surrounding the network 
is convex hull = [convex hull, convex hull [0] ] ;; the last 

vertex is the first for connections 1 sake 

area = 0 . 0D 

for i = 0, (n_elements (convexjiull) - 2) do begin 
area = area + 

( (double (network_xcors [convex_hull [i] ] ) * double (network_ycors [convex_hull [i + 
1] ] ) ) - (double (network_xcors [convex_hull [i + 1] ] ) * 
double (network_ycors [convex_hull [i] ] ) ) ) 

endfor 

area = abs (area / 2) 
areas_array = [areas_array , area] 

printf, 1, ' Area of network (in pixels A 2) : ' , 

strcompress (area) 

printf, 1, » 1 

endif 

cluster_no = cluster_no + 1 

endif 

endfor 

;; draw raster plot before it checks to see if the random distribution has 
a variance as the raster is independant of statistics 
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draw_signif icance_raster, yes_no_signif icance 
; ; draw correlations map 

//window, /free, xsize=x_size , ysize=y_size, title= ? Correl Map of All 
Many/Once Networks' ; connections__array has the data of the cells connected 

; / correl__map_image_plane, connect ions__array 
scorrel_map, symbol_array , connections_array 

; ; draw histogram of areas 

if (n_elements (areas_array) ge 2) then begin 

window, /free, title= 1 Distribution of Many/Once Areas' 
hist_areas = histogram (areas_array [1 :n_elements (areas__array) - 1], 
binsize=5000, min=0, max=max (area s_ar ray) + 5000) 

plot, lindgen(n_elements (hist_areas) + 2) * 5000L, hist_areas, 
psym=10, title= 'Distribution of Many/ Once Areas ', $ 

yrange=[0, max (hist_areas) + 1], xrange=[0, max (areas_array) + 
5000], ytitle= 'Number of Networks', xtitle='Area (Pixels A 2)' 
endif 

; ; check to see if the random distribution has any non zero values . If 
s : not, the moment cannot be defined for the distribution 
J* xmax=max (random_distribution) 

xmin=min (random_distribution) 
E3 if (xmax eq xmin) then begin ;if the min and max 

P are the same, all of the elements are the same 

ftj me s s =WIDGET_MESSAGE ( 'Random distribution has variance of zero. Try 

Sj again with greater number of iterations', /error) 
fgj print f, 1, ' ERROR i 1 

S. printf, 1, 'Moment undefined for random distribution with variance 

zero' 

close, 1 
endif else begin 

H : /plot random distribution with a line for the actual value 

0 Ip.multi = 0 

f]j bin_size = 1 

jp hist = histogram (random_distribution, binsize = 1, min = 0, max = 2 

fa * xmax + 1) /plot a histogram of the distribution 

window, /free, title='Many One Distribution' 
plot, hist, xtitle=' number of hits', ytitle= 1 frequency 1 
y2 = total (random_distribution eq true_hits) 
if y2 le 0 then begin 
/if there is no random data=true_hits , to draw the blue line 
y2 = y2 + 1 

endif 

plots, [truejiits, true_hits] , [0, y2] , color=12 /this 
draws a line where the actual number of matches lies 

/calculating the number of spikes per cell per second 
spikes_j)er_cell_per_second = (no_spikes / cells_active) / (frame_no 
* t ime__r e solution) 

/calculating the p value, standard dev. for the data 
stats = moment (random_distribution, sdev=sdev) 
no_points_right = total (random_distribut ion ge true_hits) 
p_value = no_points_right / no_iterations 

/writing more data to the file: 
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printf, 1, 'Total number of active stocks: 1 , 

strcompress (cells_active) 

printf, 1, 'Total number of frames: ', 
strcompress (frame_no) 

printf, 1, 'Total number of spikes: 
strcompress (no_spikes) 

printf, 1, 'Mean expected matches: 
strcompress (stats (0) ) 

printf, 1, 'Variance: '/ 
strcompress (stats (1) ) 

printf, 1, 'Standard deviation: 
strcompress (sdev) 

printf, 1, 'Actual matches (no. of networks): ', 
strcompress (true_hits) 

printf, 1, 'Number of networks /number of stocks: ', 
strcompress (truejiits / cells_active) 

printf, 1, 'Normalized number of networks: ', 
strcompress (truejiits / cells_active / frame_no) 

printf, 1, 'Actual/expected: 
strcompress (truejiits / stats (0)) 

printf, 1, 'Standard error for ratio: ', 
strcompress (sdev / stats (0)) 
O if (n_elements (where (size_array) ) ge 2) then begin 

Q network_stats = moment (size_array (where (size_array) ) , 

Q sdev=sdev_size_array) 

J[| printf, 1, 'Mean stocks in a network: 

\| strcompress (network_stats [0] ) 

IY| printf, 1, 'Standard deviation: ', 

7p strcompress (sdev_size_array) 
¥ * end if 

printf, 1, ' ' 

printf, 1, 'Spike firing rate: 
H : strcompress ( spikes jper_cell_per_second) 

Ql if (n_elements (distances_array) gt 2) then begin 

distance_stats = 
p moment (distances_array [1 :n_elements (distances_array) - 1], 
y , s de v= s de v_d i s t anc e s_ar r ay ) 

printf, 1, 'Mean connection distance: '/ 
strcompress (distance_stats [0] ) 

printf, 1, 'Standard deviation: ', 
s t rcompre s s ( s dev_di s t ance s_ar r ay ) 
endif else begin 

if (n_elements (distances_array) gt 1) then begin 

printf, 1, 'Connection distance: ' , 

strcompress (distances_array [1] ) 
endif 
endelse 

if (n_elements (areas_array) gt 2) then begin 

area__stats = moment (areas_array [1 :n_elements (areas_array) 
1] , sdev=sdev_areas_array) 

printf, 1, 'Mean area: ', 
strcompress (area__stats [0] ) 

printf, 1, 'Standard deviation: ', 
strcompress (sdev_areas__array) 
endif 

printf, 1, 'Significance p-value: ', 
strcompress (p_value) 
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close, 1 
endelse 

return, [random_array] 

end 



; Name : many_one_widget 

; Description: Through this widget, we can find how significant many cells 

firing in one frame is. The user 

can choose how many cells consititutes a hit and how many iterations 
the program should run through 

to create the random distribution. 



pro many_one_widget 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
f rame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 

many_one_base=WIDGET_BASE (/COLUMN, title= ' Significance of many stocks 
5 _ spiking together in one frame') 

H : many_one_f ieldl=CW_FIELD (many_one_base, /RETURNEE VENTS, / INTEGER , 

O TITLE= 'Number of stocks to spike in one frame to count as a hit:',$ 
P VALUE =4 , UVALUE=0) 

0 ge_or_eq_bgroup=CW_BGROUP (many_one_base, ['This Many or More ', 1 Exactly 

fjj, This Many'], /ROW, /EXCLUSIVE, SET_VALUE=0) 

Cj many_one_f ield2=CW_FIELD (many_one_base , / RE TURN_E VENTS , /INTEGER, 

m TITLE= 'Number of iterations : ' , $ 

P? VALUE=1000, UVALUE-0) 

T ^' ! window_size_f ield=CW_FIELD (many_one__base , /RETURN_E VENTS , /INTEGER, 

* TITLE= ' Window size for hits:', VALUE=0, UVALUE=0) 

H many_one_buttonl=WIDGET_BUTTON(many_one__base, VALUE= ' Find significance ' , 

M' UVALUE=2) 

0 WIDGET_CONTROL, /realize, many_one_base 

fjj many_one_state={many_one_f ieldl :many_one_f ieldl , $ 

many__one_f ield2 :many_one_f ield2 , $ 
y window_size_f ield: window_size_f ield, $ 

many_one_button : many_one_buttonl , $ 
ge_or_eq_bgroup : ge_or_eq_bgroup } 
WIDGET_CONTROL, WIDGET_INF0 (many_one_base , /CHILD) , 
S ET_UVALUE =many_one_s t a t e 

xmanager, 'many_one_widget ' , many_one_base 

end 



; NAME: many_one_widget_event , event 
pro many_one_widget_event , event 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
f rame^no , f i 1 e_name , $ 

total_f rame_no, time__resolution, x_size, y_size, box_size 
common with_many, filename 

stateholder3=WIDGET_INF0 (event . handler , / CHILD) 
WIDGET_CONTROL, stateholder3 , GET_UVALUE=many_one_state 
WIDGET_CONTROL, many_one__state . many_one_f ieldl , 
GET_VALUE=synchronous_cells 
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WIDGET_CONTROL , many_one__state . many_one_f ield2 , GET_VALUE==no_of_iterations 
WIDGET_CONTROL , many_one_state . window_size_f ield, GET_VALUE=window_size 
WIDGET_CONTROL , many__one_s tate . ge_or_eq__bgroup , GET_VALUE=ge_or_eq_test 

if (event. id eq many_one_state .many_one_button) then begin 

if ( (synchronous_cells le 0) OR (no_of_iterations le 1)) then begin 

mess=WIDGET_MESSAGE ( 1 Invalid fields specified!', /error) 
endif else begin 

random_array = many_cells_one_f rame (synchronous_cells , 
no_of _iterations , ge_or_eq_test , window_size) 

WIDGETMUONTROL , event . top, /hourglass 
close, 1 

openr, 1, filename 
stat_stuff = fstat(l) 
file_size = stat_stuf f . size 
close, 1 

if (file_size gt 8112) then begin 

mess = WIDGET_MESSAGE ( ' File is too large to display 
through a widget. Open it manually. If you did not save it, it is named 
Many_0ne . dat ■ , / INFORMATION) 

endif else begin 

M : xdisplayf ile, filename, title = "Statistical Data for 

Q Many/One", group = event. top, width = 75, height = 50 
0 endelse 
p endelse 
iij endif 
end 



J ! pro MultiStock 

common mother__com, pixel_array, yes_no_yalues, coef, location, cell_no, 
M ; f r ame_no , f i 1 e__name , $ 

H ! total_f rame_no, time_resolution, x_size, y_size, box_size 

Qj common old_skool_data, original__data 

fjj common flags, cells_def ined, spikes_def ined, single_plot_def ined, 
m cor re 1 coef defined 



common with_choose__cells_com, thresholdl, threshold2, rms_threshold 
common test, str_ing, state3 
common choose__cell_convert_delta, filepathl 

common textf ile_vars , text_flag, diode_array, max_num_of_diodes 



;; Init all the common block vars above! (except state3) 

pixel_array = [0] 

yes_no_values = [0] 

original_data = [0] 

coef = [0] 

location = 0 

cell_no = 1 

frame_no = 1 

f ile_name = ' ' 

total_f rame_no = 300 

time_resolution = 1.0 

x_size = 0 

y__s i z e = 0 

box_size = 0 

cells defined = 0 
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spikes_def ined = 0 
single_plot_def ined = 0 
correl_coef_def ined = 0 
rms_threshold = [0] 
thresholdl = [0] 
threshold2 = [0] 
str_ing = 'Movie. 1 
buqing_flag= 0 
diode_array = [0] 

filepathl = 'dev/null 1 ;; but where else to send our data? 
text_flag = 0 
diode_array = [0] 



Q 
P 

H 

in 

CO 



0 

III 



base=WiDGET_BASE (/column, title= 1 MultiStock ' ) 

;; setting uvalue=0 below tells the event handler that we're going to 
refer to the buttons by name, 

;; i.e. pass the button name as a string as the event. 

b=cw_bgroup (base, /row, [ ' New Experiment ' , $ 

' Save as XDR ' , $ 
'Load from XDR 1 ,$ 
1 Exit 1 ] , /return_name, UVALUE=0) 

;this is the structure that has the information for the pulldown menu 



stuff = {cw_pdmenu_s , flags :0, 

details=[ { cw_j)dmenu__s , 1, 

cw_pdmenu_s , 0 , 

cw_pdmenu_s , 2 , 

cw_pdmenu_s , 0 , 

cw_jpdmenu_s , 1 , 

cw_pdmenu_s , 0 , 

cw _pdmenu_s , 0 , 

cw_jpdmenu_s , 0 , 

cw_pdmenu_s , 0 , 

cw__pdmenu_s , 0 , 

cw_j)dmenu_s , 2 , 

cw_jpdmenu_s , 0 , 



name : ' ! } 
' Find Spikes ' } , $ 
' Root Mean Squared ' } , $ 
1 Intensity Threshold ' } , $ 
■Delete Spikes'}, $ 
'Plotting' }, $ 
'Single Plots' } ,$ 
•All Plots 1 },$ 
'All Plots-3D' } ,$ 
' Superimposed Plots ' } , $ 
'Raster Plot' } , $ 
'Overall Behavior'}, $ 
'Load from Text File'}] 



pull_down=cw_jpdmenu (base, details , /return_f ull_name , UVALUE=12) 



stuff 2= {cw_j>dmenu_s2 , flags :0 

details=[ { cw_j)dmenu_s2 , 1, 

{ cwj>dmenu_s2 , 0 , 

{ cw_j>dmenu_s2 , 0 , 

{ cw_j?dmenu_s2 , 2, 

{ cw__pdmenu_s2 , 0, 

{ cw _pdmenu_s2 , 0 , 



name : ' ' } 

'Test Significance'}, $ 
'General Significance'}, $ 
'Many Stocks One Time'}, $ 
'Two Stocks Many Times'}, $ 
'Build Correlation Map ' } , $ 
'Cross Correlogram' } ] 



pull_down2=cwjidmenu (base, details, /return_f ull_name , UVALUE=13) 



stuff3= {cw_pdmenu_s3 , flags:0, name:' 1 } 
details= [ { cwj)dmenu_s3 / 0, 'Color Tables'}, $ 
{ cw_pdmenu_s3 , 0 , ' Return to IDL ' } , $ 
{ cw__pdmenu_s3 , 0, 'Load from Excel'}] 

pull_down3=cw_j>dmenu (base , details , / return_f ull_name , UVALUE=14 ) 
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namesid=cw_f ield(base, title= ■ Filename Root '/ 
UVALUE=2, VALUE=s trying, /STRING, /return_events) 

frameid=CW_FIELD(base, title** 1 Total Number of Days 
■ , VALUE=total_frame_no, / INTEGER, / RETURN_EVENTS ) 

timeid=CW_FIELD (base, title= , Time Resolution 
' ,VALUE=time_resolution, /FLOATING, /RETURNEE VENTS) 

widget_control, /realize, base / make the widget visible 

state= { name : namesid, frame : f rameid, time : timeid} 
WIDGET_CONTROL, WIDGET_INFO (base , /Child), SET_UVALUE=state 

xmanager, 'MultiStock' , base /register this widget with the 

xmanager so it can call the event 

end 

pro MultiStock_event , event 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common old_skool_data, original_data 
*f common test, str_ing, state3 

13 common flags, eel ls_de fined, spikes_def ined, single_plot_de fined, 
0 correl_coef_def ined 

fij common with_choose_cells_com, thresholdl, threshold2 , rms_threshold 
\J common choose_cell_convert_delta, filepathl 

ff% common markov2 , yes_no_values3 , yes__no_temp , cell_no_temp , f rame_jio_temp 
5 common textfile vars, text flag, diode_array, max_num_of _diodes 

? stateholder3=WIDGET_INF0 (event . handler , /Child) 

jf* WIDGET_CONTROL, stateholder3 , GET_UVALUE==state3 ;find what is in 

^ /Child 

13 wiDGET_CONTROL, state3 . name , GE T_VALUE = f i 1 e_name 

fij WIDGET_CONTR0L, state3 . frame , GET_VALUE=total_f rame_no 

Hi WIDGET CONTROL, state3.time, GET VALUE=time resolution 

tJL. 

if ( (size (event .value) ) [1] eq 7) then begin ;; i.e. if we've pressed a 
button which returns a string as its value... 
case event. value of 

; ; FIRST ROW ITEMS 

'New Experiment 1 : begin 

end 

'Save as XDR ' : begin 

if (total_f rame_no eq 0) or (time_resolution eq 0) then 

begin 

mess=WIDGET__MESSAGE ( ' Enter frame number and time 
resolution before proceeding I 1 , /ERROR) 

endif else begin 

data__f ile=pickf ile (/write, 
f ile= ' Variables__from_exp_#_' , title= ' Create the saved variables file') 

if (data_f ile eq ' 1 ) then begin 
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mess=WIDGET_MESSAGE ( 'No saved variables file 

specified. 1 , / INFORMATION) 

endif else begin 

SAVE, /VARIABLES, FILENAME=data_f ile, all, 

/verbose 

mess=WIDGET_MESSAGE ( 'Data from this 

experiment saved! ' , /INFORMATION) 

endelse 
endelse 

end 

•Load from XDR' : begin 
state4=state3 

if (total_f rame_no eq 0) or (time_resolution eq 0) then 

begin 

mess=WIDGET_MESSAGE ( 'Enter frame number and time 
resolution before proceeding!', /ERROR) 

endif else begin 

data_file=pickf ile (/read, title= ' Select the saved 
variables file', GET_PATH=f ilepathl) 

if (data_file eq ' 1 ) then begin 
H : mess=WIDGET_MESSAGE ( 'No saved variables file 

Q specified. 1 , / INFORMATION) 

S3 endif else begin 

p rms_threshold = [0] ; ; reset it so that it 

ly can be properly reinitialized if need be 

restore, data__file 
f l size_of_loc=size (location) 

g temp_cell=size_of_loc (1) 

V* size_of_y_n_v=size (yes_no_values) 

f f rame_no=size_of _y_n_v (2 ) 

H ; cell_no=size_of_y_ji_v (1) 

if (temp_cell ne cell_no) then begin 

me s s = WIDGET_ME S SAGE ( ' Something is 
wrong with your saved variables file!', /Error) 
4 " endif else begin 

; ; cells_def ined = 1 



begin ; ; back-compatibility check 



begin 



then begin 



; ; spikes_def ines = 1 
single _jplot_def ined=0 
text_flag = 0 

if (n_elements (thresholdl) eq 1) then 

tempi = thresholdl 
temp2 = threshold2 
thresholdl = f ltarr (celljoo) 
threshold2 = f ltarr (cell_no) 
for cell « 0, cell_no - 1 do 

thresholdl [cell] = tempi 
threshold2 [cell] = temp2 

endf or 

endif 

; ; more back-compatibility. . . 

if (n_elements (rms_threshold) eq 1) 

rms_threshold = f ltarr (cell_no) 
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begin 

SET_VALUE = s t r_ing 

S ET_VALUE = f r ame_no 

SET VALUE=time resolution 



for cell = 0, cell_no - 1 do 

rms_threshold[cell] = 2.0 

endf or 

endif 

WIDGET_CONTROL, state4 .name, 
WIDGET_CONTROL , s tate4 . frame , 
WIDGET CONTROL , state4.time, 



mess=WIDGET_MESSAGE ( 'Data from 

previous experiment loaded! ' , / INFORMATION) 

endelse 
endelse 
data_f ile=0 
endelse 

end 

•Exit 1 : begin 

m WIDGET_CONTROL, /DESTROY, event . top 

m end 

n 

K ;; SECOND ROW ITEMS 

y 

m 'i 

pf 'Find Spikes. Root Mean Squared': begin 

?if; l if (total_frame_no eq 0) or (time_re solution eq 0) then 

begin 

flQ mess=WIDGET_MESSAGE ( 'Enter frame number and time 

$ resolution before proceeding!', /ERROR) 

endif else begin 

y ( if (cells_def ined ne 1) then begin 

fyl mess=WIDGET_MESSAGE ( ' You have to find cells 

before you find spikes! 1 , /error) 

endif else begin 
W rms__spikes__widget 
H ! endelse 

endelse 

end 

'Find Spikes . Intensity Threshold': begin 

if (total_jframe_no eq 0) or ( time_resolution eq 0) then 

begin 

mess=WIDGET_MESSAGE ( ' Enter frame number and time 
resolution before proceeding! 1 , /ERROR) 

endif else begin 

if (cells_def ined ne 1) then begin 

mess=WIDGET_MESSAGE( ' You have to find cells 

before you find spikes!', /error) 

endif else begin 

choose_threshold 
endelse 
endelse 

end 

'Delete Spikes': begin 

if (spikes_def ined ne 1) then begin 
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mess = WIDGET_MESSAGE ( ' You have to find spikes 
before you can delete any! 1 , /error) 

endif else begin 

delete_spikes_widget 
endelse 

end 

'Plotting. Single Plots': begin 

if (total_frame_no eq 0) or (time_re solution eq 0) then 

begin 

me ss=WIDGET_MESS AGE ( 'Enter frame number and time 
resolution before proceeding!', /ERROR) 

endif else begin 

if (spikes_def ined ne 1) then begin 
mess=WIDGET_MESSAGE ( ' You have to find spikes 
before you can plot cells!', /Error) 

endif else begin 

p ixe l_vs_ t ime_wi dge t 
endelse 
endelse 

end 

JjJ 'Plotting. All Plots': begin 

!*f if (total_f rame_no eq 0) or (time_resolution eq 0) then 

Q begin 

HI mess=WIDGET_MESSAGE ( 'Enter frame number and time 

\| resolution before proceeding! ' , /ERROR) 
fr| endif else begin 

fp if (single__plot_def ined ne 1) then begin 

mess=WIDGET_MESSAGE ( • You have to plot single 
l: cells before you can plot all cells!', /Error) 
^ endif else begin 

^ draw_spikes_many_jpages , pixel_array , 

Id yes_no_values, cell_no, frame_no, time_resolution 
fU endelse 
f!| endelse 

|s;& end 

'Plotting. All Plots-3D': begin 

if (total_frame_no eq 0) or (time_resolution eq 0) then 

begin 

mess=WIDGET_MESSAGE ( 1 Enter frame number and time 
resolution before proceeding!', /ERROR) 

endif else begin 

if (cells_def ined ne 1) then begin 

mess=WIDGET_MESSAGE ( « You have to plot single 
cells before you can plot all cells!', /Error) 

endif else begin 

draw_3D__plot , pixel_array 
endelse 
endelse 

end 

'Plotting. Superimposed Plots': begin 

if (total__frame_no eq 0) or ( time_resolution eq 0) then 

begin 
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mess=WIDGET_MESSAGE ( 1 Enter frame number and time 
resolution before proceeding !' , /ERROR) 

endif else begin 

if (cells_def ined ne 1) then begin 

mess=WIDGET_MESSAGE ( 1 You must find cells 

before you can plot them! 1 , /Error) 

endif else begin 

dr a w_sp i ke s_many_c e 1 l__de 1 1 a , p ixe l_ar r ay , 

cell_no, frame_no, time_resolution 

endelse 
endelse 

end 

1 Plotting. Raster Plot 1 : begin 

if (total_frame_no eq 0) or (time_resolution eq 0) then 

begin 

me ss=WIDGET_MESS AGE < 'Enter frame number and time 
resolution before proceeding!', /ERROR) 

endif else begin 

if (spikes_def ined ne 1) then begin 

mess=WIDGET_MESSAGE ( 1 You have to find spikes 
p before drawing a raster plot ! 1 , /Error) 
p*j endif else begin 

X if (total (yes_no_values) eq 0) then begin 

|*5 mess=WIDGET_MESSAGE ( 'None of the cells 

:^ have any spikes! Try lowering the thresholds', / INFORMATION) 
''f endif else begin 

til draw_raster 
10 endelse 
$ endelse 
|4 : endelse 
}* end 
ft 

^ ' Plotting. Overall Behavior': begin 

if (total_frame_no eq 0) or (time_resolution eq 0) then 

J"f begin 

H ; mess=WIDGET_MESSAGE (' Enter frame number and time 

resolution before proceeding!', /ERROR) 

endif else begin 

if ( single _plot__def ined ne 1) then begin 

mess=WIDGET_MESSAGE ( ' You have to plot single 
cells before you can plot all cells ! ' , /Error) 

endif else begin 

summed_sp ike s 
endelse 

endelse 
end 

'Load from Text File': begin 
end 



-THIRD ROW ITEMS - 



begin 



'Test Significance . General Significance ': begin 

if (total_f rame_no eq 0) or (time_re solution eq 0) then 
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mess = WIDGET_MESSAGE ( 1 Enter frame number and time 

resolution before proceeding I ' , /ERROR) 

endif else begin 

if (spikes_def ined ne 1) then begin 

mess = WIDGET_ME S SAGE ( 'You have to find 
spikes before calculating the significance of correlations I' , /error) 

endif else begin 

gen_sig_widget 
endelse 
endelse 

end 

'Test Significance. Many Stocks One Time':begin 

if (total_frame_no eq 0) or (time_resolution eq 0) then 

begin 

mess = WIDGET_MESSAGE ( 1 Enter frame number and time 
resolution before proceeding ! 1 , /ERROR) 

endif else begin 

if (spikes_de fined ne 1) then begin 

mess = WIDGET_MES SAGE ( 'You have to find 
spikes before calculating the significance of correlations 1 1 , /error) 
, s endif else begin 

r ' T: many one widget 

f ( % , ., — — 

endelse 

W endelse 

Q end 



'Test Significance .Two Stocks Many Times' rbegin 
p if (total_frame_jio eq 0) or (time_resolution eq 0) then 

m begin 

** mess = WIDGET_MESSAGE ( 1 Enter frame number and time 

jj\ resolution before proceeding I 1 , /ERROR) 

endif else begin 

I"* if (spikes_def ined ne 1) then begin 

(3 mess = WIDGET _MESSAGE ( 'You have to find 

rU spikes before calculating the significance of correlations 1 1 , /error) 
f| " endif else begin 

I4, two_many_widget 

endelse 
endelse 

end 

'Build Correlation Map' rbegin 

if (total_f ramejio eq 0) or (time_re solution eq 0) then 

begin 

mess = WIDGET_MESSAGE ( 1 Enter frame number and time 
resolution before proceeding !' , /ERROR) 

endif else begin 

if (spikes_def ined ne 1) then begin 

mess = WIDGET_MESSAGE ( ' You have to find 
spikes before calculating correlation coefficients!', /error) 

endif else begin 

widget__analyze 
endelse 
endelse 

end 
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'Cross Correlogram' -.begin 

if (total_f rame_no eq 0) or (time_resolution eq 0) then 

begin 

mess=WIDGET_JVIESSAGE ( 'Enter frame number and time 
resolution before proceeding I ' , /ERROR) 

endif else begin 

if (spikes_def ined ne 1) then begin 

mess=WIDGETJVIESSAGE ( 1 You have to find spikes 
before you can draw correlograms ! ' , /Error) 

endif else begin 

yes_no__values3 = yes_no__values 
choose_crosscorr 
endelse 
endelse 

end 



FOURTH ROW ITEMS 



'Color Tables ': begin 
xloadct 

end 

f% 'Return to IDL':begin 



retail 

end 

'Load from Excel' :begin 
h i 1 oad_and_c onve r t_exc e 1 f i 1 e 

CO end 



else : 
endcase 



endif 



™* end 



Ill 



: ;i f ; Object: Change pixel_vs time_widget so that you can change the threshold for 
% r *' ! certain cells 

; while not doing so for the whole array of cell. 

; This will work by having the threshold widget work as it did before, and then 
also 

; creating an option to rethreshold particular cells if necessary. 

; Through this method, the spikes will be recalculated every time the graph is 

displayed . 

; Or it can have a statement which checks whether the threshold chosen for that 
cell is the 

; universal one, and in that case it will not rethreshold. 

; WIDGET FOR DRAWING A SINGLE 

SPIKE 



pro pixel_vs_time_widget 

common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common with_choose_cells__com, thresholdl, threshold2, rms_threshold 

deep_base=WIDGET_BASE (/row, title='Plot of a Single Stock', /scroll) 
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lef t_base=WIDGET_BASE (deep_base , /column) 
right_base=WIDGET_BASE (deep_base , / column) 
draw=WIDGET_DRAW(left_base, xsize=450, ysize=450) 
;; kludge to fix the 'only one cell selected' bug 
if (cell_no gt 1) then begin 

slide=WIDGET_SLIDER(left_base, value=l, maximum^celljio, minimum=l / $ 

title= ' Choose the stock you want to plot', 

uvalue=' slider_pressed' ) 
endif else begin 

slide=WIDGET_TEXT (lef t__base, value='Only one stock selected...') 

endelse 

spike_max=CW_FIELD(right_base, title='Y Max for Spike :' , VALUE=0, 
/ RETURN_E VENTS , /FLOATING) 

spike_min=CW_FIELD(right_base, title='Y Min for Spike: \ VALUE=-20, 
/ RE TURN_JE VENT S , /FLOATING) 

ymax=CW_FIELD (right_base, title='Y Range Max:', 
VALUE=max( [transpose (pixel_array (0 , *) ) , 10]), / RE TURN_E VENTS , /FLOATING) 

ymin=CW_FIELD(right_base, title='Y Range Min: 
VALUE=min( [transpose (pixel_array(0, *) ) , -10]), /RETURN_E VENTS, /FLOATING) 

xmax=CW_FIELD (right_base, title='X Range Max:', 
VALUE=total_frame_no*time_resolution / / RE TURN_E VENTS , /FLOATING) 
u xmin=CWJ?IELD(right_base, tit le= » X Range Min : ' , VALUE=0, / RETURNEE VENT S , 

IZ /FLOATING) 

W median_filter_window=CW_FIELD (right _base, title= ' Window for Median Filter: 

0 » , VALUE=2 0, / RE TURN_E VENTS , /FLOATING) 

0 smoothing_window==CW_FIELD(right_base, title= 1 Window for Mean Smoothing:', 

III VALUE = 3 , /RETURNEE VENTS, /FLOATING) 

\j buttonsl=cw__bgroup (right _base, /row, ['Smooth', 'Median Filter ',' Restore 

(ft Original Waveform'] , $ 

BUTTON_UVALUE= [ 1 smooth_jpressed ' , ' f ilter_pressed ' , 
^ 'restore_pressed' ] , UVALUE=5) 

f rms threshold box=CW FIELD(right base, title='RMS Threshold:', VALUE=2.0, 

!*$ " 

^ ! UVALUE=19, / RE TURN_E VENTS , /FLOATING) 

0 rms_spikes_type_bgroup=CW_BGROUP (right_base, ['Positive Spikes', 'Negative 

111 Spikes'], /ROW, /EXCLUSIVE, SET_VALUE=0) 

0 ' threshold_one=CW_FIELD{right_base, title='2 Frame Intensity Diff 
|4 : Threshold: \ VALUE =thres ho ldl [0] , $ 

UVALUE=6, /RETURN_E VENTS, /FLOATING) 
threshold_two=CW_FIELD (right Joase, title='3 Frame Intensity Diff 
Threshold: ' , VALUE=threshold2 [0] , $ 

UVALUE=7 , / RE TURN__E VENTS , / FLOATING) 
buttons2=cw_bgroup (right_base, /row, ['Plot with RMS', 'Plot with 
Intensity Difference'], $ 

BUTTONJJVALUE= [ ' RMS_plot_pressed ' , 1 dif f _j>lot_pressed ' ] , UVALUE=3) 
delete_spike=CW_FIELD(right_base,title=' Delete Spike Number:', VALUE= ' 0 ' , 
UVALUE=2 , $ 

/ RETURNEE VENTS , / INTEGER) 
buttons 3 =cw_bgroup (right_base, /row, ['Delete Stock', 'Delete All Trailer 
Spikes ' ] , $ 

BUT TON_UVALUE = [ ' delete_cell jpressed ' , 
' delete_trailer_spikesjpressed' ] , UVALUE=4) 

holder^ { draw : draw, s lide : s lide , delete_spike : delete_spike , $ 
spike_max: spike_max, spikejnin: spike_min, $ 
ymax:ymax, ymin:ymin, xmin:xmin, xmax:xmax,$ 
rms_spikes_type_bgroup : rms_spikes_type_bgroup , $ 
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median_f ilter_window:median_f ilter_window, 

smoothing_window: smooth ing_window, $ 

rms_threshold_box : rms_thresholdjDOX, threshold_one : threshold_one , 

threshold_two : threshold_two} 

WIDGET_CONTROL , de ep_bas e , /realize 
WIDGET_CONTROL, deep_base, set_uvalue=holder 

plot, findgen(frame_no) * time__re solution, pixel_array (0 , *) , 
yrange=[min( [transpose (pixel_array (0 , *) ) , -10]), max ( [transpose (pixel_array ( 0 , 

*>>, 10])],$ . n 

xrange=[0, total_f rame_no * time__resolution] , title=l, ytitle='% 

Change in Index', xtitle= ' Time (days)', FONT=-l 

for frames=l, frame_no - 1 do begin ;this doesn't 

consider a spike at the 1st frame 

if (yes_no_values (0, frames) ne 0) then begin 

plots, [frames * time_re solution, frames * time_resolution] , 

[0, -20] 

endif 

endf or 

XMMTAGER, 1 pixel_ys_time_widget ' , deep_base 

end 

n 

O 

0 " 

fil pro pixel__vs_time_widget_event , event 

\| common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
Ip f rame_no , f i le_name , $ 

}fj total_frame_no, time_resolution, x_size, y_size, box_size 

common share_w_j)ixel_vs_time_widget_com, y_min, y_max, spike_min, spike_max, 
jf, x_max, x_min 

common with_choose_cells__com, thresholdl, threshold2, rms_threshold 
^; common old_skool_data, original__data 

P common flags, cells__def ined, spikes_def ined, single_plot_de fined, 
fU correl_coef_def ined 

|4 ; WIDGET_CONTROL , event . top , get_uvalue=state 

WIDGET_CONTROL , event . id , GET JUVALUE=uval 
WIDGET_CONTROL , state. draw, GET_VALUE= window 
WIDGET_CONTROL, state . slide, GET_VALUE=number_l 

if (string (uval) eq ' slider_pressed 1 ) then begin 

WIDGET_CONTROL, state . delete_spike , set_value = 0 
spike_l = 0 

WIDGET_CONTROL, state. ymax, set_value = max (pixel_array ( (number_l - 



1), *)) 
1), *)) 



y_max = max (pixel_array ( (number_l - 1), *) ) 

WIDGET_C0NTR0L, state. ymin, set_value = min (pixel_array ( (number_l - 
y_min = min (pixel_array ( (number_l - 1) , *) ) 
; ; check for 10,-10 boundary. . . cheap. . . 

WIDGET_CONTROL , state. ymax, set_value = max([y_max, 10]) 
y_max = max([y_max, 10]) 

WIDGET_CONTROL, state. ymin, set_value = min([y_min, -10]) 
y_min = min([y_min, -10]) 
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WIDGET__CONTROL , state . rms_threshold_box, 
SET_VALUE=rms_threshold [number_l - 1] 

rms_t = rms_threshold[nuTnber_l - 1] 

WIDGET_CONTROL, state . threshold_one , SET_VALUE=thresholdl [number_l - 

1] 

tl = thresholdl [number_l - 1] 

WIDGET_CONTROL, state . threshold_two , SET_VALUE=threshold2 [nuntoer_l - 



1] 



endif 



t2 = threshold2 [number_l 



1] 



WIDGET_CONTROL , state . delete_spike , GET_VALUE=spike_l 

WIDGET JTONTROL, state .median_f ilter_window, GET_VALUE=median_f ilter_window 

WIDGET_CONTROL, state . smoothing_window, GET_VALUE=smoothing_window 

WIDGET_CONTROL, state . rms_threshold_box, GET_VALUE=rms_t 

WIDGET_CONTROL , state . threshold_one , GET_VALUE = 1 1 

WIDGET~CONTROL, state . threshold_two, GET_VALUE=t2 

WIDGET_CONTROL, state .ymin, GET_VALUE=y_min 

WIDGET_CONTROL , state . ymax, GET_VALUE=y_max 

WIDGET_CONTROL , state. xmin, GET_VALUE=x_mi n 

WIDGET_CONTROL , state.xmax, GET_VALUE=x_max 

WIDGET_CONTROL, state . spike_min, GET_VALUE=spike_min 

WIDGET_CONTROL, state . spike_max, GET_VALUE=spike_max 



si 

m 



wset, window 

windows, this will still plot in the original window 



;even if you open other 



m 
id 



; ; trailing spike clumps deletion 

if (string (event. value) eq ' delete J:railer_spikes_pressed' ) then begin 
for frame_counter = (frame_no - 1), 2, -1 do begin 

if ( (yes_no__values [number_l - 1, f rame_counter] eq 1) and 
( (yes_no_values [number J - 1, f rame_counter - 1] eq 1) or 
(yes_no_values [number_l - 1, f rame_counter - 2] eq 1) ) ) then begin 

yesjio_values [numberJL - 1, f rame_counter] = 0 

endif 

endfor 

draw_spikes, number__l 

endif 

;/ if we're re-plotting, update the yes_no_values to current threshold 
if (string (event. value) eq 1 dif f _j>lot_pressed' ) then begin 

WIDGET_CONTROL , state . delete_spike , set_value = 0 

spike_l = 0 

if ((tl eq 0) OR (t2 eq 0)) then begin 

mess=WIDGET_MESSAGE ( ' Threshold cannot be 0!', /Error) 

WIDGET_CONTROL , state . threshold_one , 
SET_VALUE=thresholdl [number_l - 1] 

WIDGET_CONTROL , state . threshold_two , 
SET_VALUE=threshold2 [number_l - 1] 

tl = thresholdl [number_l - 1] 

t2 = threshold2 [number_l - 1] 
endif else begin 

thresholdl [number_l - 1] = tl 

threshold2 [number_l - 1] = t2 

yes_no_values=make_single_binary (pixel__array , yes_no_values , 
number_l, frame_no, thresholdl [number_l - 1], threshold2 [number_l - 1]) 
endelse 
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endif 



if (string (event .value) eq 1 RMS_plot_jpressed ' ) then begin 
WIDGET_CONTROL, state . delete_spike , set_value = 0 
spike_l = 0 

WIDGET_CONTROL , state . rms_spikes_type_bgroup , GET_VALUE=pos_or_neg 
rms ^threshold [number_l - 1] = rms_t ;; update new threshold with 
value in textbox 

if (pos_or__neg eq 0) then begin 

cell_stats = moment (pixel_array [number_l - 1,*], 

s de v= eel l_df f _s t dde v ) 

cell_df f_mean = cell_stats [0] 

for frame_counter = l r frame_no - 1 do begin 

if (pixel_array [number_l - 1, f rame_counter] gt 
(cell_dff_mean + (rms_threshold [numberJL - 1] * cell_df f_stddev) ) ) then begin 

yes_no_values [number_l - 1, f rame_counter] = 1 
endif else begin 

yes_no_values [number_l - 1, f rame_counter] = 0 
endelse 

endf or 

I j ; ; negative spikes is pos_or_neg = 1 

IZ endif else begin 

W cell_stats = moment (pixel jar ray [number_l - 1,*], 

^ ! sdev=cell_df f_stddev) 

Q cell_dff_mean = cell_stats [0] 

Jjj for frame__counter = 1, framejao - 1 do begin 

\j if (pixel_array [number_l - 1 , f rame_counter] It 

m (cell_df f_mean - (rms__threshold [number_l - 1] * cell_df f_stddev) ) ) then begin 
p5| ~" yes_no__values [number_l - 1, f rame^counter] = 1 

endif else begin 

?. yes_no_values [number_l - 1, f rame_counter] = 0 

endelse 

^ endfor 

•per, 

M endelse 
III endif 

y ; ;; filtering- subtractive median filter: see f ilter_medxan_subtractive .pro 

if (string (event. value) eq ' f ilter_pressed ' ) then begin 

if ( (median_f ilter__window It 2) or (median_jf ilter_window ge 
f rame_no) ) then begin ;; error checking 

error_message = WIDGET_MESSAGE ( 1 Window size for filtering must 
be greater than 1 and less than the number of f rames ! ' , /INFORMATION) 
endif else begin 

f ilter_median_subtractive, pixel_array [number_l - 1,*], 
median__f ilter_window, f iltered_cell 

pixel_array [numberJL - 1,*] = f iltered_cell 
endelse 

endif 

; ; smoothing, via the IDL function 'smooth' 

if (string (event. value) eq 1 smooth_pressed' ) then begin 

if ( (smoothing_window It 2) or ( smoothing_window ge f rame_no) ) then 
begin ;; error checking 

error_message = WIDGET_MES SAGE ( 1 Window size for smoothing must 
be greater than 1 and less than the number of f rames !' , /INFORMATION) 
endif else begin 



NY02:354955.1 



-78- 



pixel_array [numberJL - 1,*] = smooth (pixel_array [number JL 
1 , * ] , smoothing_window) 
endelse 

endif 

; ; back to basics 

if (string (event .value) eq ' restore ^pressed' ) then begin 

pixel_array [number_l - 1, *] = original_data [number_l - 1, *] 



1) , *)) 
l), *)) 



endif 



WIDGET_CONTROL, state . ymax, set_value = max (pixel_array ( (number_l - 
yjnax = max (pixel__array ( (number_l - 1); *) ) 

WIDGET_CONTROL, state. ymin, set_value = min (pixel_array ( (number_l - 
y_min = min (pixel_array ( (number_l - 1), *) ) 
;; check for 10,-10 boundary... cheap... 

WIDGET_C0NTR0L , state. ymax, set_value = max([y_max, 10]) 
y__max = max ( [y__max, 10] ) 

WIDGET_CONTROL , state. ymin, set_value = min([y_min, -10]) 
y_min = min([y_min, -10]) 



Ml] 

10 

63 1] 



J:;:f ;; correct for zero-threshold anyways, it ! s cheap 

O if ( (tl eq 0) OR (t2 eq 0)) then begin 

mess=WIDGETjyiESSAGE ( 1 Threshold cannot be 0! 1 , /Error) 
WIDGET_CONTROL , state . threshold_one , SET_VALUE=thresholdl [number_l 

WIDGET_CONTROL, state . threshold_t wo , SET_VALUE=threshold2 [number_l 

endif 

; ; update the list of spikes 

time_array = where (ye s_no_values (number_l - 1, *) eq 1) ;this finds 

^ the location of the spikes for each cell 
I'll number spikes = n_elements ( time_array) 

\>& ; ; check for spike deletion 

if ( (time_array (0) ne -1) AND (spike_l ne 0)) then begin 

if ((spike_l gt number_spikes) OR (spike_l It 1)) then begin 

mess=WIDGET_MESSAGE ( 1 You entered an invalid spike number. 1 ) 
WIDGET J20NTROL, state . delete_spike , set_yalue = 0 
spike_l = 0 
endif else begin 

yes_no_values (number_l - 1, time_array (spike_l -1)) =0 
draw_spikes, number_l ; that cell 

will now be displayed 

WIDGET J^ONTROL, state . delete_spike , set_value = 0 

spike_l = 0 
endelse 

; ; here is where all the plotting gets done 

endif else begin ;if no spikes are to be 

deleted 

draw_spikes , number_l 

WIDGET_CONTROL , state . delete_spike , set_value = 0 
spike_l = 0 
endelse 
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; ; check for cell deletion 

if (string (event. value) eq 1 delete_cell_jpressed ' ) then begin 
yes_no_values (number__l - 1,*) = 0 
draw_spikes , number_l 

mess=WIDGET_MESSAGE (string ( ' Stock ' , byte (number_l) , ' has been 
successfully deleted. ') , /information) 
endif 

single_plot_defined=l ;this is sent to Multicell_event 

so the other programs now that single plots has been done 

e nd . 

; Name: rms_spikes_widget 

; Description: a new, improved method of spike detection! Finds the mean and a 
number of standard deviations 

(g enera lly 2... the Root Mean Square version of signal/noise handling), 

and any points beyond this 

threshold are considered spikes. 



ji^pro rms_spikes_widget 

common rms spikes, rms__spikes_state 

|* rms_spikesjDase=WIDGET_BASE (/COLUMN, title='Find Spikes- RMS') 

W rms_spikes_button=WIDGET_BUTTON(rms_spikesjDase, VALUE= ' Find Spikes', 

*^UVALUE=2) 

' h 4 rms_spikes_type_bgroup=CW_BGROUP (rms_spikes_base , [ ' Positive 

©Spikes' , 'Negative Spikes'], /ROW, /EXCLUSIVE, SET_VALUE=0) 

|T| rms_spikes _threshold_f ield=CW_FIELD (rms_spikes_base , / RE TURN__E VENT S , 

i; / FLOATING, TITLE= 'Threshold (Number of Std Devs) ',$ 
U VALUE=2.0, UVALUE=0) 

S'7 rms__spikes_median_filter_window_size_field-CW_FIELD(rms_spikes_base, 
IZ /RETURN_E VENTS, /FLOATING, TITLE= ' Window Size for Median Filter ',$ 
p VALUE=2 0.0, UVALUE=0) 

rms spikes text=WIDGET TEXT (rms spikes base, VALUE= ' Set window size to 0 
to prevent filtering. . . ' ) 

U 

WIDGET_CONTROL, /realize, rms_spikes_base 

rms_spikes_state= { rms_spikes_base : rms_sp ike Sybase , $ 

rms_spikesjDutton:rms_spikes_button, $ 
rms_spikes_type_bgroup : rms_spikes_type_bgroup , $ 

rms spikes_threshold__f ield: rms__spikes_threshold_f ield, $ 

rms_spikes_median_filter_window_size_field:rms_spikes_median_filter_window 
_size__f ield} 

WIDGET_CONTROL, WIDGET_INFO (rms_spikes_base , /CHILD) , 
SET_UVALUE = rms_sp i ke s_s t a t e 

xmanager , 1 rms_spikes_widget ' , rms_spikes_base 

end 

; Name: rms_spikes_widget_event 

pro rms_spikes_widget_event , event 
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common mother_com, pixel_array; yes_no_values , coef, location, cell^no, 
frame_no, file_name,$ 

total_frame_no, time_resolution, x_size, y_size, box__size 
common flags, cells_def ined, spikes_de fined, single_plot_def ined, 
correl_coef_def ined 

common with_choose_cells_com, thresholds threshold2, rms_threshold 
common rms_spikes, rms_spikes_state 

rms J3tate_stuf f =WIDGET_INFO (event . handler, /Child) 

WIDGET_CONTROL, rms_state_stuf f , GETJJVALUE=rms_j3pikes_state ;find 
what is in /Child 

if (event. id eq rms_spikes_state . rms_spikes_button) then begin 
;; read interface data... 

WIDGET_CONTROL , rms_spikes_s tate . rms_spikes_type_bgroup , 
GET__VALUE=pos_or_neg 

WIDGET_CONTROL, rms_spikes_state . rms_spikes_threshold_f ield, 

GET_VALUE=threshold 

WIDGET_CONTROL, 

rms_spikes_state . rms_spikes_median_f ilter_window_size_f ield, 

s , GET VALUE=window size 

p Si - - 

M for i = 0, cell_no - 1 do begin 

0 rms_threshold[i] = threshold 

t! endfor 

ru 

%i if (window_size gt 0) then begin 

ffj f ilter_median_subtractive, pixel_array, window_size, 

fp;pixel_array 

end 

yes no values=intarr (cell_no, f rame_no) 

h& ~ ~ 

P ;; positive spikes is pos_or_neg = 0 

111 if (pos_or_neg eq 0) then begin 

p " for cell_counter = 0, cell_no - 1 do begin 

l2 cell_stats = moment (pixel_array [cell_counter , *] , 

sdev=cell_df f_stddev) 

cell_dff__mean = cell_stats [0] 

for frame_counter = 1, frame_no - 1 do begin 

if {pixel__array [cell_counter , f rame_counter] gt 
(cell_df f_mean + (threshold * cell_df f_stddev) ) ) then begin 

yes_no_values [cell_counter, f rame_counter] = 

1 

endif 

endfor 

endfor 

; ; negative spikes is pos__or_neg = 1 
endif else begin 

for cell_counter = 0, cell_no - 1 do begin 

cell_stats = moment (pixel_array [cell_counter, *] , 

s dev=c e 1 l_df f _s t ddev ) 

cell_df f_mean = cell_stats [0] 

for f rame__counter= 1, frame_no - 1 do begin 

if (pixel_array [cell_counter , f rame__counter] It 
(cell_df f_mean - (threshold * cell_df f_stddev) ) ) then begin 
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yes_no_values [cell_counter , f rame_counter] 

endif 



endf or 

endf or 
endelse 



endif 

end 



spikes_de fined = 1 

mess=WIDGET_MESSAGE( 'Spikes have been found!', / INFORMATION) 
WIDGET_CONTROL, rms_spikes_state . rms_spikes_base, /DESTROY 



; NAME: 

SCORREL_MAP 
; PURPOSE: 

This procedure creates a correlation map between the input cells. A 
correlation i, 

a circular representation of correlation between cells of the slice. Those 
.. cells that 

r A : are correlated in either direction are joined by lines which are 

^proportional to their 

Q; correlation coefficients 

p; PARAMETERS: 

fjj; symbols: this is the array of names of companies. 

\j; output or coef_array: this is the correlation coeff. array 

plpro scorreljnap, symbols, output 

* window, /free, xsize=500, ysize=500, title= 1 Correlation Map Between Stocks* 

N' 

H i number=n_elements (symbols) 

P ?a DRAWS CIRCLE OF THE 

f|J CELLS 

m radius = 5* (number A 2) /through trial and error, this seems 

1^ the best radius 

of f setl=f ltarr (number+1) 



of f set2=f ltarr (number+1) 
of f set3=f ltarr (number+1) 
of f set4=f ltarr (number+1) 

plot, f ltarr (10), $ /need to call this blank plot procedure 

just to set the ranges for the xyouts 

xrange= [-8*number^2 , 8 * number A 2 ] , $ 

yrange= [-8*number A 2, 8*number A 2] , $ 

xstyle=4, ystyle=4, xmargin= [0, 0] , ymargin= [0, 0] 

erase /this is so you 

don't see any junk left over.. 

FOR 1=0, NUMBER- 1 DO BEGIN 

ANGLE = ( ( !pi*2) /number) *i 

of f setl (i) =cos (angle) 

of f set2 (i) =sin (angle) 

of f set3 (i) = (radius*1.25) *cos (angle) 

of f set4 (i) = (radius*l .25) *sin (angle) 

xyouts, l.l*radius*offsetl (i) , 1 . l*radius*of f set2 (i) , symbols (i) , 

/ alignment 
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endf or 



; PRINT 

LINE S 

? output_array=print_chi (cut) 

; this part prints the lines between highly correlated pairs 
-this is done since the output we have has an empty cell at postion 7. 
rthis finds the cell numbers from thier positions 
for s=0, number- 1 do begin 

for t = 0, number-1 do begin 

;if (output (s, t) gt 0) AND (output (s,t) ne 1) then begin 
if (output (s,t) gt 0) then begin 
plots, [radius*of f setl (s) , radius*of f setl (t) ] , 
[radius*of f set2 (s) , radius*of f set2 (t) ] , thick=3*output (s , t) , linestyle=2 

; arrow, off setl (a), offset2(s), offsetl(t), offset2(t), /data, 
thick=3 *array ( s , t ) 

endif 

endf or 

endfor 
H; !p.font=0 
Cj end 

u 

o 

flj pro summed_spikes 

\j common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
Ipi f r ame_jio , f i 1 e_name , $ 

l£ total frame no, time resolution, x size, y size, box size 

%Q - - - - - - 

J base=WIDGET_BASE (/column, title= 1 Overall Stock Behavior 1 ) 

sigma=CW_FIELD (base, title= f Enter the value for sigma (smoothing 
M factor) : ' , VALUE=5, $ 

P UVALUE = 1 , / RETURN_E VENT S , / FLOAT ING ) 

flj buttonl=widget__button (base , value= f PLOT 1 , uvalue= ' plotjressed 1 ) 

f| WIDGET_C0NTR0L, base, /realize 

| ;; |, widget_control , base, set_uvalue= sigma 

xmanager, 1 summed_spikes 1 , base 

end 



pro summed_spikes_event , event 

common mother_com, pixel_array, yes_no_values , coef, location, cell__no, 
f rame_no , f i le_name , $ 

total_f rame_no, time_re solution, x_size, y_size, box_size 

W IDGET__CONTROL , event . top , get_uvalue= sigma 
WIDGET_C0NTR0L , s igma , GET_VALUE= s igma 
window, 8 , title= ' Stock Behavior' 

1 except =0 

summed_array=total (yes_no_values, 1) 
h=dblarr (f rame_no) 
j=dblarr (frame_no) 
z=dblarr (frame no) 
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for t=0., (frame_no -1), 1 do begin 
h[t]=0 

j [t]=o 

for s=0, (framejio -1) do begin 

for i=l,Summed_array [s] do begin ; tne 

two below divisions are unnecessary » . , . , , 

h[t]=h[t] + (exp(l))"(((--5 * (s - t) A 2) ) / sigma 2) / 

(Sqrt(2 * !pi * sigma A 2)) 
endfor 

endf or 

for i=0, (frame_no-l) do begin 

j[t]=j[t] + (exp(l)) A (((-.5 * (i - t) A 2) ) / sigma 2) / 

(Sqrt(2 * ipi * sigma A 2) ) 
endfor 

z[t]=h[t]/j [t] 

endfor 

plot / findgen(frame_no)*time_resolution /Z/ title= 'Plot of overall stock 
i :! ;.behavior • , $ 

^ xtitle = 'Time (days)', ytitle = 'h' 

!except=l 

#end 

0 



J y ; Name : two_cells_many_times . pro 

Nl; Description: This procedure looks at the cells that fire together more than 
CO once. It then counts the number of times that pairs 

;6; of cells have fired more than once together. It also does this for random 
/cells that have been created through the procedure 

\„k) RANDOM_TEST . The random cells are tested for multiple hits through the 



y procedure MULTIPLEJTESTJSIGNIFICANCE. 
m l} Explanation of Variables: no_iterations is the number of iterations the 
^program will loop to create the random distribution. 

least no of matches is the fewest number of times two 



1: 



J cells have to spike together to begin counting the correlation. 

H window_size is how far, to left and to right, of a spike 

we look for other spikes for two cells to be considered coactive 

pro two_cells_many_times, least_no_of ^matches , no_iterations , window_size 

common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
frame_no, file_name,$ 

total_frame_no, time__resolution, x_size, y_size, box_size 
common with_create_dist , cells_active ;this contains the 

number of active cells 

forward_f unction multiple_test_signif icance 
common with_two, filename 

common stockdata, symbol_array , date_array 

connections^array = bytarr (cell_no , cell_no) 
yes_no_signif icance = yes_no_values 

/initializing the counters to zero 
true hits = 0 
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/creating the random distribution 
no_spikes = total (ye s_no_values) 

random_distribution = intarr (no_iterations) ;this array 

holds one random number of matches for every iteration 

for t = 0, no_iterations - 1 do begin ; to 

create the random dist. we have to repeat the 2 steps no_iterations times 

make_random_data, seed, r andom_ar r ay , num_active_cells 

random_distribution(t) = count_random_hits_2_manyX (random_array , 

least no of matches, window_size, num_active_cells) 

— — j. — /makes 
endfor 

an array with a number for each iteration 
number corresponds to the number of times 2 cells 

least 'least no of_times ' in the random movie for the t-th iteration 



the 

fire at 



/printing some of the data 

rilename= pickfile (/write, f ile= ( ' Two_Stocks__' + 
strcompress (string (leas t_no_of_matches) ) + 'JTimes' ) ) 
if (filename eq ' ' ) then begin 

mess=WIDGET_MESSAGE('This data will not be saved! /INFORMATION) 
f ilename= 1 Two __Many . dat 1 

end if 
close, 1 

openw, 1, filename 
printf, 1, 'Statistical Data:-' 
11 printf, 1, 'Number of times two stocks must spike together to count as a 

Njhit : ' , least_no_of ^matches 

ffl printf, 1, 'Number of iterations: ', no__iterations 

printf, 1, 'Total number of active stocks: « , cells_active 



m 

rte' 



*. /finding the actual number of matches 

f7 distances_array = [0.0] 

^ for cell_l_counter = 0,cell_no - 2 do begin 

13 for cell_2__counter = cell_l_counter + 1, cell_no - 1 do begin 

flJ tempjiits = 0 

temp_locations_cell_l = [-1] 
temp_locations_cell_2 = [-1] 

for window_counter = -window_size , window__size do begin 
temp_cell_l = intarr (frame_no + (2 * window_size) ) 
temp_cell_2 = intarr (frame_no + (2 * window_size) ) 
temp_cell__l (window_size : f rame_no+window_size-l) = 
yes_jio_values (cell_l_counter , *) 

temp_cell_2 (window_size : f rame_no+window_size-l) = 
yes_no_values (cell_2_counter , *) 

temp_cell = (temp_cell_l * shif t (temp_cell_2 , 

window_counter) ) 

; ; ;temp_cell = (yes_no_values (cell_l_counter, *) * 
shift (yes_no_values (cell_2_counter , *) , window_counter) ) 

if ( (size (where (temp_cell) )) [0] eql) then begin 
temp_hits = (temp_hits + total (temp_cell) ) 
temp_locations_cell_l = 
[temp_locations_cell_l, where (temp_cell) - window_size] 

temp_locations_cell_2 = 
[temp_locations_cell_2, where (temp_cell) - window_counter - window__size] 

endif 

endfor 
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if (temp_hits ge least_no_of_matches) then begin 
;; rearrange location arrays: 

;; ...we want to first get rid of all -1 elements 
(meaning no paired spikes) , then sort the array, then (as sorting returns array 
subscripts 

;; and not actual elements) we need to determine what 
the original values for the new subscripts are. If you're still confused, 

;; run this line of code in sections on test data, e.g. 
'a = [-1,1,2,-1,4,8,-1,5] & a (sort (a (where (a ge 0) ) ) ) « . 

real_locations__cell_l = 
(temp_locations_cell_l (where (temp__locations_cell_l ge 
0) ) ) (sort (temp_locations_cell_l (where (tempJLocations_cell_l ge 0) ) ) ) 

real_locations_cell_2 = 
(temp_locations_cell_2 (where (temp_locations_cell_2 ge 

0))) (sort (temp_locations_cell_2 (where ( temp_locations_cell_2 ge 0)))) 

; ; compute connection distance 
connection_distance = 
sqrt ( ( (double (location (cell_l_counter) . coord [0] ) - 
double (location (cell_2_counter) .coord [0] ) ) 2) + 
, . ( (double (location (cell_l_counter) . coord [1] ) - 
double (location (cell_2_counter) .coord [1] ) ) 2) ) 



CI distances_array = [distances_array , connect ion_di stance] 

13 connections_array (cell_l_counter, cell_2_counter) = 1 

Q truejiits = true_hits + 1 

j'jj ; / dump to file 

printf, 1, 'Stock number one and location: 



ffj strcompress (cell_l_counter + 1), location (cell_l_counter) . coord 
lZ printf, 1, 'Frames cell one spikes in:' 

^ printf, 1, real_locations_cell_l 

J printf, 1, 'Stock number two and location: 

^ : strcompress (cell_2_counter + 1), location (cell__2_counter) . coord 



H printf, 1, 'Frames cell two spikes in:' 

0 printf, 1, real_locations_cell_2 

111 printf, 1, ' ' 

fl ;; update 'significance array' which is used to draw the 



|4 raster plot 

for f rame_counter = 0, temp_hits - 1 do begin 
ye s_no_s i gni f i c anc e ( c e 1 l__l_c ount e r , 
real_locations_cell_l [f rame_counter] ) = TEMPORARY (yes_no_signif icance (cell_l_count 
er, real_locations_cell_l [f rame_counter] ) ) + 1 

yes_no_signif icance (cell_2_counter, 
real_locations_cell_2 [f rame_counter] ) =TEMPORARY (ye s_no_signif icance (cell_2_count 
er, real_locations_cell_2 [f rame_counter] ) ) + 1 

endf or 

endif 

endf or 

endf or 

;draw raster plot before it checks to see if the random distribution has a 
variance as the raster is independant of statistics 
draw_signif icance_raster, ye s_no_si gni f icance 

; check to see if event occurs randomly and thus if a p value can be 
calculated. 

xmax=max (random_distribution) 
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xmin==min (random distribution) 



if (xmax eq xmin) then begin 

mess=WIDGET_MESSAGE ( 1 Random distribution has variance of zero. Try 
again with greater number of iterations', /error) 
printf, 1, 1 ERROR I ' 

printf, 1, 'Moment undefined for random distribution with variance 

zero ' 

free_lun, 1 
close, 1 
p_value=-l 
endif else begin 

;making the histogram for the random distribution 

!p.multi=0 

bin_size=l 

hist=histogram (randomjiistribution, binsize=l , min=0 , 
max= (2*xmax) +1) ;plot a histogram of the distribution 

window, /free, title='Two Many Distribution' 
plot, hist, xtitle=' number of hits', ytitle= ' frequency ' 
y2=total (random_distribution eq true_hits) ; f or drawing the 

blue line for the real data, we have to find the height of that line in the 
^histogram 

Q if y2 le 0 then begin ;if there are no random 

pvalues equal to the true data, draw a line of height one 

5 y 2=1 

ft endif 

• : 1 plots, [true_hits, truejiits] , [0, y2] , color=12 ;this draws a 

5 .'?line where the actual number of matches lies 

sill 



m 



/calculating the number of spikes per cell per second 
spikes _j>er__cell_per_second= (no_spikes/cells_active) / (f rame_no*time_resolut 



3:f. rlOn 

h 

"AW 

m 

u 
u 



/draw the correlation map 

connections_array= ( connect ions_array) +TRAJSTSPOSE (connections_array) 
/this is to make the connections_array symmetric. Because if cell 1 is 
connected to 2 , 2 is also connected to 1. correl_map_plane needs the symmetric 
array 

scorrel_map, symbol_array , connect ions_array 

//window, /free, xsize=x_size, ysize=y_size, title= 1 Correl Map for 



Two Many' 



/ / correl_map_image_plane , connections_array 



/find the p value, standard deviation etc. 
stats-moment (random_distribution, sdev=sdev) 
noj?oints_right=total (randomjiistribution ge true_hits) 
p_value=no_points__right/no__iterations 



printf, 1, 
printf, 1, 
printf, 1, 
strcompress (stats (1) ) 
printf, 1, 



/print the data to a file 
printf, 1, 'Total number of frames: 
Total number of spikes : 
Mean expected matches : 
Variance : 



'Standard deviation: 



strcompress (f rame_no) 
strcompress (no_spikes) 
strcompress (stats (0) ) 
i 

1 , strcompress (sdev) 
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printf, 1, 'Actual matches: ' , strcompress (truejiits) 

printf, 1, 'Actual/expected: ' , 

strcompress (true_hits/ stats (0) ) 

printf, 1, 'Standard error for ratio: 
strcompress (sdev/stats (0) ) 

printf, 1, 'Spike firing rate: f , 
strcompress (spikes_per_cell_per_second) 

if (n_elements (distances_array) gt 2) then begin 
distance_stats = 
moment (distances_array [1 :n_elements (distances_array) -1] , 
sdev=sdev_distances_array) 

printf, 1, 'Mean connection distance: 1 , 

strcompress (distance_stats [0] ) 

printf, 1, 'Standard deviation: 
strcompress (sdev_distances_array) 
endif else begin 

if (n_elements <distances_array) gt 1) then begin 

printf, 1, 'Connection distance: ', 
strcompress (distances_array [1] ) 
endif 
endelse 

printf, 1, 'Significance p-value: ! , strcompress (p_value) 

close, 1 
endelse 



b 

s- 

fl ;i - 

; Name : two_many__widget 

N ; Description: Using this widget the user can test the significance of two 

W cells firing many times in 

|0 ; the movie to analyze. The user will have to input the 

i specifications such as the number of times 

jj.^ ; that two cells fire together to be taken as the min. criterion for a 

| ; | : hit. 

rjgpro two_many__widget 

common mother_com, pixel_array, yes_no_values, coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_f rame__no, time_re solution, x_size, y_size, box_size 

two_many_base=WIDGET_BASE (/ COLUMN, title= ' Significance of two stocks 
spiking together many times') 

fieldl=CW_FIELD(two_many_base, / RE TURN_E VENTS , / INTEGER, TITLE='Two stocks 
should spike together at least this many times:', VALUE=2 , UVALUE=0) 

field2=CW_FIELD(two_many_base, / RETURNJE VENT S , /INTEGER, TITLE= 1 Number of 
iterations:', VALUE=1000, UVALUE=0) 

field3=CW_FIELD(two_many_base, / RETURNEE VENT S , /INTEGER, TITLE= ' Window 
size for hits:', VALUE =0 , UVALUE=0) 

buttonl=WIDGET_BUTTON(twojnany__base, VALUE= ' Find significance', UVALUE=2) 

WIDGET_CONTROL, /realize, two_many_base 

two_many_state={f ieldl : f ieldl , f ield2 : f ield2 , f ield3 : f ield3 , 
but tonl : buttonl } 

WIDGET_CONTROL, WIDGET_INFO ( two_many_base , /CHILD) , 
SETJJVALUE=two_many_state 

xmanager , ' two_many_widget ' , two_many_base 

end 
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NAME : two_many_widget_event 



pro two_many_widget_event , event 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_frame_no, time_re solution, x_size, y_size, box_size 
common with_two, filename 

two_many_state_inf o=WIDGET_INFO (event . handler , /CHILD) 
WIDGET_C0NTR0L, two_many_state_inf o , GETJCJVALUE=two_many_state 
WIDGET_CONTROL, two_many_state . f ieldl, GET_VALUE=least_no_of_matches 
WIDGET_CONTROL, two_many_state . f ield2 , GET_VALUE=no_of_iterations 
WIDGETJTONTROL, two_many_state . f ield3 , GET_VALUE==window_size 

if (event. id eq two_many_s t ate .butt onl) then begin 

if ( (least_no__of jnatches le 0) or (no_of_iterations le 1) or 
(window_size It 0) or (window_size gt f rame_no) ) then begin 

mess=WIDGET__MESSAGE (' Invalid fields specified!' , /ERROR) 

endif else begin 

two_cells_many_times , least_no_of _matches , no_of _iterat ions , 



window_size 



WlDGET_CONTROL , event . top , /hourglass 
close, 1 

% openr, 1, filename 

5;;f stat_stuf f =f stat (1) 

W f ile_size=stat_stuf f . size 

fU close, 1 

% 4 if (file_size gt 8112) then begin 

flj mess=WIDGET_MESSAGE ( 'File is too large to display 

tQj through a widget. Open it manually. If you did not save it, it is named 
I " Two JVIany . da t 1 , / INFORMATION) 

endif else begin 

f7 xdisplayfile, filename, title = "Statistical Data for 

P>Two Many", group = event. top, width = 75, height = 50 
JJ* endelse 



$ endelse 

% 

H : end 



W endif 



; Widget for creating correlation coefficients matrix and drawing general 
correlation maps. 

pro widget_analyze 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 

f r ame_no , f i 1 e_name , $ 

total_f rame^no, time_resolution, x_size, y_size, box_size 
basel=WIDGET_BASE (/COLUMN, title= ' General Slice Correlation') 
type= { cw_pdmenu_s , flags : 0 , name : 1 1 } 

nitty_gritty =[ { cw_j)dmenu_s , 0, 'Calculate Correlation Coeff Matrix'}, 

$ 

{ cwj)dmenu_s, 1, 'Build Correlation Map'}, $ 
{ cw_jpdmenu_s , 0 , 'On Image ' } , $ 
{ cw_pdmenu_s, 2, 'Spatial'}] 
pull_down2=CW_PDMENU (basel, nitty_gritty , /RETURN_FULL_NAME) 

WIDGET_CONTROL, basel, /realize 
xmanager, ' widget_analyze ' , basel 
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end 



; Name: widget_analyze_event 

; Description: This procedure is called when you press the 1 Analyze' button on 
the main menu. 

; It allows you to calculate the correlation coefficient (WHICH MUST BE THE 
FIRST STEP) and then 

; draw a correlation map on an image or spatial . 
pro widget_analyze_event , event 

common mother_com, pixel_array, yes_no_values, coef, location, cell_no, 
f r ame_no , f i 1 e_name , $ 

total_frame_no, time__resolution, x_size, y_size, box_size 

common flags, cells_def ined, spikes_def ined, single _j>lot_def ined, 

correl_coef_def ined 

CASE (event .value) of 

'Calculate Correlation Coeff Matrix': BEGIN 
choose_correl 

END 

|J 'Build Correlation Map. On Image': BEGIN 

W if (correl_coef_def ined ne 1) then begin 

y mess=WIDGETjyiESSAGE ( ' You have to first calculate the 

^correlation coefficients! ' , /Error) 

>h 4 endif else begin 

fg scorrel_map 

endelse 

END 

'V- 

J* 'Build Correlation Map . Spatial ' : BEGIN 

y if (correl_coef_def ined ne 1) then begin 

fO mess=WIDGET_MESSAGE ( ' You have to first calculate the 

P correlation coefficients! 1 , /Error) 
Ms endif else begin 

bas e 5 =WIDGET_BASE (/COLUMN, title= ' Spatial Correlation Map') 

draw7=WIDGET_DRAW (baseS , xsize=x_size , ysize=y__size) 

WIDGET_CONTROL , /realize, baseS 

correl_map_image_j?lane , coef 

endelse 
END 

ENDCASE 
end 
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APPENDIX B 



; to find argument of the maximum value of the array 
; used by hidden_markov_event .pro 

function arg_max, x,y,j,t,D,A 

Z=dblarr (y+1) 

for i=x,y do begin 

if ( A[i,j] eq 0) then begin 

z [i] =-10D307 

endif else begin 

z[i] = D[t-l,i] +alog(A[i, j] ) 

endelse 

endf or 

m=max (Z , k) 
return, k 

end 

; function used by HMM to find the probability of observation vector of 
y,; length cell_no, at time t, given the array of Poisson lambdas B, and statejao 

^function B_prob, state_no, observation, t, B, cell_no 

;f ; 1 except =2 

ftJs=double(l) 

if; J 

fflfor i=0, cell_no -1 do begin 
f fli ob s e rva t i on ( i , t ) = doub le{observation(i,t)) 
J^if (B [state_no, i] eq 0) then begin 
L y.if (observation (i,t) eq 0) then begin 

l^endif else begin 
*fj|s=s * 0 
l5f endelse 

W endif else begin 

|^s=s* ( (exp (-B [state_no, i] ) ) * ( (B [state_no, i] ) A (observation (i , t) ) ) / 
factorial (observation (i , t) ) ) 

;s=s - B[state_no,i] + (observation (i, t) ) *alog (B [state_no, i] ) - 

alog (factorial (observation (i , t) ) ) 

endelse 

endf or 

; s=exp (s) 

return, s 

end 

pro choose_init _J?ar 



common mother_com, pixel_array, yes_no_values , coef , location, celljio, 
frarae__no, file__name, direct, $ 

total_f rame_no, time_re solution, x__size, y_size, box_size 
common markov, 

n, m, A, B , P , observation, f , g, d, Fi , Q, back, xsi , Gamma, n_i , a_i , b_i , p_i , c , w_s , t_max 
common indexes , i_a , i_b , i_j> 
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i_a=0 
i_b=0 
ijp=0 

Values^ [ ' random' , 'uniform'] 

Values2= [ 1 random' , 1 uniform ',' averaged over intervals'] 
base=widget_base (/column, title= 1 Initial parameters generation 1 ) 

Ainit=widget_droplist (base, value=values, uvalue=0, title= ' Choose initial state 
transition prob. A: ') 

Binit=widget_droplist (base, value=values2 , uvalue=l, title= ' Choose initial state 
characteristcs B: ') 

Pinit=widget_droplist (base, value=values , uvalue=2 , title= ' Choose initial state 
prob . P : ' ) 

Apply=widget_button (base , value= • Apply 1 , uvalue=3 ) 

widget_control, base, /realize 
xmanager, ' choose_init_par ' , base 

end ___ 

pro choose_init_par_event , event 



^common mother_com, pixel_array, yes_no_values, coef, location, cell_no, 
p£f rame_no, file_name, direct, $ 

!j?total_f rame_no, time_resolution, x_size, y_size, box_size 

j;j ( fcommon markov, 

I'Un, m, A, B , P , observation, f , g, d, Fi , Q , back, xsi , Gamma, n_i , a_i , b_i , p_i , c , w_s , t_max 
^common indexes, i__a,i_b,i_p 

CO 

fyiWIDGET_CONTROL, event, id, GET_UVALUE=uval 



||!!lF (TAG NAMES (event, / STRUCTURE NAME ) EQ ' WIDGET_DROPLIST ' ) $ 

Q 



THEN BEGIN 
Hi Case uval of 
0 0 : begin 

M= CASE event . index OF 



0 : begin 

; print, 'A random' 
i_a=i_a+l 

; Random distribution 

for i = 0, (n-1) do begin 

for j=0, (n-1) do begin 

A [i, j ] =randomu (seed) 

endf or 

endf or 

zl=dblarr (n) 

zl=total (A, 2) 

for i=0, (n-1) do begin 

for j=0, (n-1) do begin 

A[i, j]=A[i, j]/zl[i] 

endf or 

endf or 

end 
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1 : begin 

; print, ' A uniform 1 
i_a=i_a+l 

/Uniform distribution 

for i=0, (n-1) do begin 

for j=0, (n-1) do begin 

A[i, j]=l/n 

endf or 

endf or 

end 



END CASE 

end 
1 : begin 

CASE event. index OF 



0 : begin 

/print, 1 B random' 
i_b=i_b+l 

for 1=0, (n-1) do begin 
for j=0, (m-1) do begin 
B [i, j ] =Randomu (seed) /Random 
endf or 
endf or 



zl= total (b,2) 

for i=0, (n-1) do begin 

% 4 for j=0, (m-1) do begin 

CO B[i,j]=B[i,j] / zl[i] ; think 

Cp endfor 

$ endfor 

jU; 

m 1: begin 

*p* /print, ? B uniform' 

^ ib=i b+1 

ffy — — 

W ; uniform 

M< for i=0, (n-1) do begin 

for j=0, (m-1) do begin 

; mm=m . 

B[i,j]=l/float(m) 

endfor 

endfor 

end 



2 : begin ; average over uniform segments 
i_b=i_b+l 

for j=0, m-1 do begin 
for i=0, n-1 do begin 
s=0 
s_i=0 

for k=f ix(i* ( (frame_no - w_s + l)/n)) , f ix ( (i+1) * ( (f rame_no -w_s + 
l)/n)-l) do begin /change frame_no 
s=s+observation [ j , k] 
s_i=s_i+l 
/print, k 
endfor 

/print, 'dupa' 
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s=double (s) 

s_i=double (s__i) 

B [i , j ] =double ( s / s__i) 

endf or 

endfor 

end 

END CASE 



end 



2 : begin 

CASE event . index OF 

0 : begin 

; print, 'P random' 
i_p=i_p+l 

for i=0, (n-1) do begin 
P ( i ) =randomu (seed) ; random 
endfor 
z2=total (p) 

for i=0, (n-1) do begin 
J* P(i)=P(i)/z2 
C3 endfor 
%rS end. 

III 1 : begin 



\$ ; print, 'P uniform' 

CO 

II] 



IB i_p=i_p+l 
; uniform 
P(0)=1 

; 5: for i = l, (n-1) do begin 

M« P(i)=0 

H» endfor 

0 end 

f[j END CASE 



end 
endcase 
ENDIF 



if (uval eq 3) then begin 

print, 'Apply pressed' 

; print, 1 i_a= ' , i_a, ' , ijb= f # i_b, ! , i_p 

if (i_a eq 0) then begin /make random 
; print, 'A random' 
for i=0, (n-1) do begin 
for j=0, (n-1) do begin 
A[i, j] =randomu (seed) 
endfor 
endfor 
zl=dblarr (n) 
zl=total (A, 2) 
for i=0, (n-1) do begin 
for j=0, (n-1) do begin 
A[i, j]=A[i,j]/zl[i] ;?? 
endfor 



NY02:354955.1 



endf or 

endif 

if (i_b eq 0) then begin ; make random 
/print, 'B random' 
for i=0, (n-1) do begin 
for j=0, (m-1) do begin 
B [ i , j ] =Randomu (seed) 
endf or 
endf or 

zl=total (b,2) 

for i=0, (n-1) do begin 

for j=0, (m-1) do begin 

B[i, j]=B[i, j] /zl[i] /think 

endf or 

endf or 

endif 

if (ij> eq 0) then begin / make random 
/print, ! P random' 
for i=0, (n-1) do begin 
P (i) =randomu (seed) 
endfor 
z2=total (p) 

for i=0, (n-1) do begin 
P(i)=P(i)/z2 
endfor 
Jendif 

Jjprint, 'These are initial state transition probabilities, A' 
Sprint, A 
fnprint , 1 space ' 

Sprint, 'These are initial probabilities of observing symbol m at state n, B 1 

? print , b 

f " ! print, ' space ' 

j^print, 'This is the initial state distribution, P ! 
imprint, p 
ftjprint, ' space 1 

M ; n_i=n /for save function 
A_i=A 
B_i=B 
P_i=P 

widget_control, event .top, /destroy 
endif 

end 

pro hidden_markov 



□ 



common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
frame_no, file_name, direct, $ 

total_f rame_no, time_resolution, x_size, y__size, box_size 
common markov, 

n,m, A,B, P, observation, f ,g,d, Fi , Q, back, xsi , Gamma, n_i , a__i , b_i , p_i , c , w_s , t_max 
common markov2 , yes_no_values3 , yes_no_temp , cell_no_temp , f rame_no_temp 



/ this is the definition of a test input with 400 frames and 4 states 
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; f rame_no=40 0 

;observation=intarr (f rame_no) /observation sequence 

;for i=0,99 do begin 

;x=randomn (seed, binomial= [1, .3] ) 

;if (x eq 1) then begin 

/observation [i] =1 

;endif else begin 

/observation [i] =2 

/ endelse 

/ endf or 

/for i=100,199 do begin 
;x=randomn(seed,binomial= [1, .6] ) 
/if (x eq 1) then begin 
/observation [i] =1 
/endif else begin 
/observation [i] =3 
/endelse 
/ endf or 

/for 1=200,299 do begin 

/x=randomn(seed,binomial= [1, . 1] ) 

/if (x eq 1) then begin 
„ . / observation [i] =2 
J/ endif else begin 
^1; observation [i] =0 
0; endelse 
CI /endf or 

flJ/ for i=300,399 do begin 
M/x=randomn (seed, binomial= [1, .5] ) 

if (x eq 1) then begin 
|ji;y=randomn(seed,binomial= [1, .4] ) 
*'*;if (y eq 1) then begin 
f a /observation [i] =1 
f*} endif else begin 
f*; observation [i] =2 
C|/ endelse 
||j/ endif else begin 

fl;y=randomn (seed, binomial= [1, .4] ) 
iy./if (y eq 1) then begin 

/observation [i] =0 

/endif else begin 

/observation [i] =3 

/ endelse 

/ endelse 

/ endf or 

/m=double (4) 



/this transforms yes_ no_values to observation seq., information is lost 

/because only one cell per time unit is allowed to spike 

; observation=intarr (f rame_no) /observation sequence 

/for i=0, (frame_no -1) do begin 

/observation [i] =0 

/for j=0, (cell_no -1) do begin 

/if (yes_no_values [ j , i] eq 1) then begin 

/observation [i] =j+l /needs change 

/ endif 
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; endf or 
; endf or 

;m=double (cell_no + 1) ; number of observation symbols (+1 if neither cell fires 
- 0) 



;this transforms yes_no_values to observation by adding everything - creates 
artificially many zeros 



3J f ;; 

a 



; m=double (cell_no + 1) ; number of observation symbols (+1 if neither cell fires 
- 0) 

observationl=intarr (f rame_no * cell_no) /observation sequence 
z=0 

for j=0, (frame_no -1) do begin 
for i=0, (cell__no -1) do begin 
if (yes_no_values [i, j ] eq 1) then begin 
observationl [z] =i+l 
endif else begin 
observationl [z] =0 
endelse 
z=z+l 
endf or 
endf or 

frame no=frame no*cell no 



m 

JU;this creates observation from yes-no_values by adding everything but only if 

• ^cells fire simulatneously 

* ; m=double (cell_no + 1) ; number of observation symbols (+1 if neither cell fires 

H/observation_2=intarr (f rame_no * cell_no) /observation sequence 
z = 0 

; ; for j=0, (frame_no -1) do begin 
z2 = 0 

for i=0, (cell_no -1) do begin 
if (yes_no_values [i, j ] eq 1) then begin 
observation^ [z] =i+l 
z=z+l 
z2 = l 

endif else begin 
observation_2 [z] =0 
endelse 
endf or 

if (z2 eq 0) then begin 
z-z+1 
endif 
endf or 

; observation2=intarr (z) 
;for i=0, (z-1) do begin 
;observation2 [i] =observation_2 [i] 
; endf or 
; f rame__no=z 
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; observat ion=observat ion2 



test=fltarr (2,3) 
test [0,0] =1.1 
test [0,1] =1.2 
test [0,2] =1.3 
test [1,0] =2.1 
test [1,1] =2.2 
test [1,2] =2.3 
suml=f ltarr (5) 
suml=total (test, 1) 
sum2=f ltarr (5) 
sum2=total (test, 2) 
print, test 
print, 'this is sum 1' 
print, suml 
print, 'this is sum 2' 
print, sum2 

|1jbase=widgetJoase (/column, title= 'Hidden Markov Modelling' 
\Jstuff= { cw_pdmenu_s , flags :0, name: ' 1 } 
;!fjdetails= [{ cw_pdmenu_s , 1, 'Create input data' }, $ 

^ t a _ r\ i 

CO 
u 



cw__pdmenu_s , 0 , 

cw_pdmenu_s , 0 , 

cw_j?dmenu_s , 0 , 

cw_j?dmenu_s , 2 , 

cw__pdmenu_s , 0 , 

cw _pdmenu_s , 0 , 



},$ 



1 Save init. pars . ' , 
/return name, 



from observation sequence '}, $ 

'from shuffled in time observation sequence '}, $ 
'from Monte Carlo simulated observation seq (space) 
'test' } , $ 

1 Record y__n_values 1 } , $ 
'Retrieve y_n_values ' } ] 
f : i|pull_down=cw_j)dmenu(base, details, /return_name, UVALUE=4) 
pbuttonl=cw_bgroup(base, /row, ['Choose initial parameters.', 
' Load init. pars . 1 , 'Save final pars.', 'Load final pars . ' ] , 
UVALUE=0) 

state_no=CW_FIELD (base, title= 'Enter the number of states you would like to use 
in the model :» , VALUE=4 , $ 

UVALUE=5 , / RETURNEE VENTS , /Floating) 
but ton2=widget_but ton (base, value= ' Estimate model parameters.', uvalue=3) 
iteration_no=CW_FIELD (base, title=' Enter the number of iterations you would like 
to perform:', VALUE=2 0, $ 

UVALUE=4 , / RETURN_E VENT S , /FLOATING) 
button3=widget_but ton (base, value='Find the most probable hidden state 
sequence. ' , uvalue=2) 

button4=widget_button(base, value= ' Compute P (observation | model parameters) 1 , 
uvalue=l) 

button5=widget_button(base, value='Find cross correlations within states', 
uvalue=6) 

correlation_state=CW_FIELD (base, title='In which state do you want to find 
cross-correlations ?:', VALUE=0, $ 

UVALUE=7, /RETURNEE VENTS , /FLOATING) 
state={ state__no : state_no, iterationjno : iterationjno, 
correlation_state : correlation_state} 
widget_control , base, set_uvalue=state 
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WIDGET_CONTROL, statejno, GET_VALUE =n 
widget_control , base, /realize 

xmanager, ' hidden_markov' , base 
end 



pro hidden_markov_event , event 

common mother_com, pixel_array, yes_no_values , coef , location, cell_no, 
frame_no, file_name, direct, $ 

total_frame_no, time_resolution, x_size, y_size, box_size 
common markov, 

n, m, A, B, P, observation, f ,g,d,Fi ,Q, back, xsi, Gamma, njL, a_i,b_i,p_i, c, w_s, t_max 
common markov2 , yes_no_yalues3 , yes_no_temp, cell_no_temp, f rame_no_temp 

; !except=2 
n_old=n 

WIDGET_CONTROL , event . id , GET_UVALUE=uval 

Widget_control, event .top, get_uvalue=state 

WIDGET__C0NTR0L , state . iteration_no, GET_VALUE=iteration_jao 

^ W I DGET_ CONTROL , state . Stat e_no , GET_VALUE=n 

pjw I DGET_ CONTROL , state . correlationjstate , GET_VALUE=correlation_state 

m-. 

WCASE uval OF 

rii 

H|4 : Begin 



frlCase event. value OF 

from observation sequence' : BEGIN 



I is 

^ /print , 1 these are y_n_values 1 
M ; print , yes_no_values 



fU ; print , 



? ; m-cell_no 

w_s= 30. ; size of sliding window, step is equal to 1 

observation=intarr (cell_no, FRAME_N0 - w_s +1) 

for i=0, (frame_no - w_s ) do begin 

for k=0, (cell_no -1) do begin 

s=0 

if (i It (frame_no - w__s +1)) then begin 
for j=i, (i+w_s -1) do begin 
s=s+yes_no_values [k, j ] 
endf or 

observation [k, i] =s 

endif 

endfor 

endf or 

print, 'Input from observation sequence created' 
end 



'from shuffled in time observation sequence': BEGIN 
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yes_no_values2=intarr (cell_no, f rame_no) 
shuf f le=intarr (cell_no) 

for 1=0, cell_no -1 do begin 

shuf f le [i] = (f rame_no / 2) * randomn (seed, uniform=l) 
endf or 

for i=0, cell_no -1 do begin 

for j=0, frame_no -1 - shuf fie [i] do begin 

yes_no_values2 [i, j] =yes_no_values [i, j+shuffle [i] ] 

endfor 

for j=(frame_no - shuf fie [i] ) , (frame_no -1) do begin 
yes_no_values2 [i, j ] =yes_no_values [i, shuf fie [i] - frame_no + j] 
endfor 
endfor 

/print, 'These are yes and no values: 1 
; print , ye s_no_value s 2 
L| ; ; print, ' ' 

Hi ; ye s_no_va 1 ue s =y e s__no_va lue s 2 
pon=cell_no 

lr w — s= 30 • ; size of sliding window, step is equal to 1 
?^observation= intarr (cell_no, FRAME_NO - w_s +1) 
Nfor i=0, (frame_no - w_s ) do begin 
Cflfor k=0, (cell_no -1) do begin 

ias=o 

c if(i It (frame_no - w_s +1)) then begin 
| ( |,for j=i, (i+w__s -1) do begin 
y.s=s+yes_jio_values [k, j ] 
J^. endf or 

^observation [k, i] =s 
jjjfendif 
)4 endfor 
I 4 ' endfor 

print, 'Input from shuffled in time observation sequence created' 
end 



'from Monte Carlo simulated observation seq (space)': BEGIN ;this is copied 
from vikram's make__random_data .pro 



;find the number of spikes each cell of the true data has 
no_spikes=intarr (cell_no) /this is the array with the number of 

spikes for each cell 

for cell=0, cell_no-l do begin 

no_spikes (cell) =total (yes_no_values (cell, *) ) 

endfor 

cells_active=total (no_spikes gt 0) 

/create a random binary spike train. First, numbers for locations of 
spikes are generated randomnly making 

/sure that no number is repeated more than once. Next, a binary spike 
train is created using these locations 
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/for the position of spikes. Each cell has the same number of spikes as 
the true data had. 

random_array=intarr (cell_no, f rame__no) 
for cell=0, cell_no-l do begin 
index=no_spikes (cell) -1 

; using a random number generator for a uniform distribution 
;here we see if there are any spikes in the original cell of the 

data 

if index ge 0 then begin 
repeat begin 

temp=randomu (my_seed, no_spikes (cell) ) ; random 

numbers 

temp=f ix (f rame_no*temp) 

;this has the location of the unique elements in the 

random row of data 

;this is the test to make sure no numbers are repeated 

in the data 

unique= (uniq (temp (sort (temp) )) ) ;# of 

y /unique elements 

jj>*& numberl =n__elements (unique) /total number 

Sof elements 

*J number 2 =n_elements (temp) 

W endrep until numberl eq number2 

flj 

S3 random_ar ray (cell, temp-l)=l ;put spikes at 

fHlthe random locations 

endfor 

| :?! i i yes_no_values=random_array 
I; /print , 1 this are y__n values 1 
p /print, y e s_no_va lue s 

Hfm=cell_no 

C3w_s= 30. / size of sliding window, step is equal to 1 
H'observation=intarr (cell_no, FRAME_N0 - w_s +1) 

for i = 0, (frame_no - w_s ) do begin 

for k=0, (cell_no -1) do begin 

s=0 

if (i It (frame_no - w_s +1)) then begin 
for j=i, (i+w_s -1) do begin 
s=s+yes_no_values [k, j ] 
endfor 

observation [k, i] =s 

endif 

endfor 

endfor 

print, 'Input from monte carlo simulated data created' 
end 

'test ' : BEGIN 

; this is going to be a test 

print, 'Frames 1-100 should be in state 0 with paramters : [ 3,2,8,1]' 
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print, 'Frames 101-200 should be in state 1 with parameters: [ 3,4,8,4] ' 
print, 'Frames 201-300 should be in state 2 with parameters: [ 3,3,6,0] ' 
print, 'Frames 301-400 should be in state 3 with parameters: [11,3,6,4]' 



m=4 . 

cell_no=4 . 

w_s=ll 

f rame_no=410 

observation=intarr (cell_no, 
for i=0, (99 ) do begin 
observation 
observation 
observation 
observation 
endfor 
for i=100, 
observation 
observation 
. observation 
observation 



FRAME NO - w s +1) 



0, i] =randomn(seed, poisson=3) 

1, i] =randomn (seed, poisson=2) 

2 , i] =randomn (seed, poisson=8) 

3,i] =randomn (seed, poisson=l) 

199 ) do begin 

0, i] =randomn (seed, poisson=3) 

1, i] =randomn (seed, poisson=4) 

2 , i] =randomn (seed, poisson=8) 

3, i] =randomn(seed, poisson=4) 

2 99 ) do begin 

0 , i] =randomn (seed, poisson=3) 

1, i] =randomn (seed, poisson=3) 

2,i] =randomn(seed, poisson=6) 
3,i] =0 



;;;; endfor 
Of or i=200, 
Ms observation 
|1| observation 
\| observation 
| : p ob s e rva t i on 

endfor 
J**for i=300, 
f 3 observation 

observation 
^ ob s e rva t i on 
fl-3 observation 
III endfor 
P print, 'test input created 

3 ; 

end 



3 99 ) do begin 

0,i] =randomn (seed, poisson=ll) 

1, i] =randomn(seed, poisson=3) 

2 , i] =randomn (seed, poisson=6) 

3 , i] =randomn (seed, poisson=4) 



'Record y_n_values': begin 

yes_no_temp=yes_no__yalues 

eel l_no_temp=cell_no 

f r ame_no_t emp = f r ame_no 

print, 'Yes and no values recorded' 

end 



'Retrieve y_n_values ' : begin 

cell_no=cell_no_temp 

f r ame_no = f rame_no_t emp 

yes_no_values=intarr (cell_no, f rame_no) 

ye s_no_va lue s =y e s_no_t emp 

print, 'Yes and no values retrieved' 

end 

' test2 ' : begin 
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for i=0,4 do begin 
print, yes_no_values [i, 7] 
endfor 
end 



else : 
endcase 

/print, 'This is the input:' 
; print , observation 
/print, ' ' 

t_max=f rame_no - w_s +1 
n=double (n) 
A=dblarr (n,n) 
B=dblarr (n,m) 
P=dblarr(n) 



end 



0 : Begin 
s:Case event. value OF 



Choose initial parameters.': BEGIN 

Q 

C3 
fij 

\Jif (n_old ne n) then begin 
|fljn=double (n) 
£QA=dblarr (n,n) 
B=dblarr (n,m) 



y P=dblarr (n) 
f' ; endif 

choose_jLnit jar 

HI end 

0 

H= 'Save init. pars.': BEGIN 

data_file=pickf ile (/write, f ile= ' HM_IVars_exp ' , title= ' Create the 
saved variables file') 

if (data_file eq ' ' ) then begin 

mess=WIDGET_MESSAGE ( 'No saved variables file specified' , 

/ INFORMATION) 

endif else begin 

Save, FileName=data_f ile, A_i, B__i, P_i, n_i 
mess=wlDGET_MESSAGE (' Initial data from this experiment saved 1 , 

/INFORMATION) 

endelse 



end 



'Load init. pars.': BEGIN 

data__f ile=pickf ile (/read, title= ' Select the saved 
variables file', GET__PATH=f ilepathl) 

if (data_file eq 1 ' ) then begin 
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mess=WIDGET_MESSAGE ( 1 No saved variables file specified', 

/ INFORMATION) 

endif else begin 
restore, data_file 

mess=WIDGET_MESSAGE ( 1 Previous initial A, B, P, and n 

variables loaded' , / INFORMATION) 
endelse 

data_f ile=0 

a=a_i 

b=b_i 

P=P_i 
n=n_i 

print, 'This is initial A' 

print, A 

print, ' ' 

print, 'This is initial B 1 

print, B 

print, 1 1 

print, 'this is initial p' 

print, p 

print, 1 1 



WIDGET_CONTR0L, state . state_no, SET_VALUE=n 



0 

enc j 

K 6 

% y 

V| "Save final pars.': BEGIN 

m 

|| data_file=pickfile (/write, f ile- ' HM_FVars_exp ' , title- 1 Create the 

*/ saved variables file') 

i, if (data_file eq * ' ) then begin > _ 

f; mess=WIDGET_MESSAGE ( 'No saved variables file specified 1 , 

J* / INFORMATION) 

W endif else begin 

flj Save, FileName=data_f ile, A, B, P, n 

0 mess=WIDGET_MESSAGE ( 'Data from this experiment saved', 

\£ /INFORMATION) 

endelse 

end 

'Load final pars.': BEGIN 

data_file=pickf ile (/read, title= ' Select the saved 
variables file', GET_PATH=f ilepathl ) 

if (data__file eq ' * ) then begin 
mess=WIDGET_MESSAGE ( 'No saved variables file specified', 



/INFORMATION) 



/ INFORMATION) 



endif else begin 
restore , data_f ile 

mess=WIDGET_MESSAGE ( ' A, B, P, and n variables loaded', 



endelse 
data_f ile=0 

print, 'This is A' 
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print, A 
print, 1 1 

print, 'This is B 1 
print, B 
print, 1 ' 

print, 'this is p' 
print, p 
print, ' ' 

WIDGET_CONTROL, state . state_no , S E T_VALUE = n 

end 



else : 

endcase 

end 

1: Begin 

I , ;The Log probability of the observation given the model 

^ g=double (0) 

;;;;f ;for i=0, n-1 do begin 

;g=g+f (t_max-l, i) 
U ; endf or 

RJ for t=0,t_max -1 do begin 
Nl 9=9 ~ aio9 (c (t) ) 

fpi endf or , , n T 

print, 'This is the Log probability of the observation sequence gxven the model 

^ print, g 
: print, 1 1 
s . end 

flj ;2 : begin ; Viterbi with no scaling 
Q ; print, 'button 2 pressed 1 

; This is an implementation of the viterbi algorithm, the pupose of 

; viterbi algorithm is to find the single most probable state sequence given 

; the observation sequence and model parameters. 

; iexcept=2 

;D=dblarr (t_max,n) /highest probability along a single state path of the 
observation sequence 

;Fi=dblarr(t_max,n) ; the backtracking step to retrieve the hidden states 
; initiation 

;for i=0, n-1 do begin 

;D [0, i] =p [i] *B_prob (i , observation, 0) 

;Fi [0, i] =0 

; endf or 

/recursion 
; Z=dblarr (n) 

/for t=l, t_max-l do begin 
/for j-0, n-1 do begin 
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for i-0, n-1 do begin 
Z[i]=D[t-l,i]*A[i,j] 
endf or 

D [t, j] = (max(Z) ) *B [ j , (observation (t) ) ] 

Fi[t, j]==arg_max(0,n-l, j,t,D,A) ; i changes from 0 to n-1 
endf or 
endf or 

termination 
for i = 0, n-1 do begin 
Z [i] =D [t_max-l, i] 
endf or 

Prob_path=max(Z,k) ; probability of the given path 
Q=intarr (t_max) ; Q is the hidden state sequence 
Q[t_max -l]=k 

for t=t_max-2 ,0,-1 do begin 
Q[t]=Fi[t+l, (Q[t+l])l 
endf or 

print, 'This is D 1 
print, D 

print, 'this is Fi ' 
print, Fi 
print, 'This is Q' 
print, Q 
print, ' ' 

, print, 'This is the corresponding probability of the hidden state sequence 
Itj! given the observation data:' 
print, Prob_path 
print, 1 ' 
end 



fi 



'■Hi 

CO 



> u 2 : begin 

if* ; print, 'button 2 pressed' 

0 ; This is an implementation of the viterbi algorithm, the pupose of 

W ; viterbi algorithm is to find the single most probable state sequence given ^ 

|3 • the observation sequence and model parameters, here i use logarithmic scaling 

D=dblarr (t_max,n) ; highest probability along a single state path of the 
observation sequence 

Fi=dblarr (t_max,n) ; the backtracking step to retrieve the hidden states 

/ initiation 

for i=0, n-1 do begin 

if ( (P(i) eq 0) or (B_prob (i, observation, 0 , B, eel l_no) eq 0) ) then begin 

D [0, i] =-10D307 

Fi [0, i] =0 

endif else begin 

D [0, i] =alog (P (i) ) +alog (B_jprob (i , observation, 0 , B, cell_no) ) 
Fi [0, i] =0 
endelse 
endf or 



; recursion 
Z=dblarr (n) 

for t=l, t_max-l do begin 
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for j=0, n-1 do begin 

for i=0, n-1 do begin 

if ( A[i,j] eq 0) then begin 

Z [i] =-10D307 

endif else begin 

Z [i] =D[t-l,i] +alog (A[i, j] ) 

endelse 

endfor 

if ( B_jprob(j, observation, t,B,cell_no) eq 0) then begin 
D[t, j]=-10D307 

Fi[t,j]=arg_max(0,n-1, j,t,D,A) ; i changes from 0 to n-1 
endif else begin 

D[t, j] = (max(Z) ) +alog(B_prob(j observation, t,B,cell_no) ) 

Fi[t,j]=arg_max{0,n-1, j,t,D,A) ; i changes from 0 to n-1 

endelse 

endfor 

endfor 

; termination 
for i=0, n-1 do begin 
Z [i] =D [t_max-l, i] 
U endfor 

^ ln_of JProb_path=max(Z,k) ; probability of the given path 

Q=intarr (tjnax) ; Q is the hidden state sequence 
d. Q[t_max -1] =k 
O for t=t_max-2, 0, -1 do begin 
til Q[t]=Fi[t+l, (Q[t+13)] 
h} 4 endfor 

|P print, 'This is Q* 
III print, Q 

- print, 1 ' , 

lU print, 'This is LogP of the hidden state sequence given the observation data: 

l r : print, ln_of_Prob_jpath 

J£ print, 1 ' 

^ end 

Q 3 : begin 

H= print, 'button 3 pressed' 

for i-0, iteration_no do begin 

par__reestimation 

endfor 

print, 'this is new A' 
print, A 
print, ' ' 

print, 'This is new B 1 
print, B 
print, ' ' 

print, 'this is new P 1 
print, p 
print, 



end 

6 : begin 
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; cell_no=2 
;W_S=3 

; f rame_no=15 

/ t_max=frame_no - w_s + 1 

;yes_no_values=intarr (cell_jio, f rame_no) 

/yes_no_values [0 , *] =f indgen (f ratne_no) 

;yes_no_values [1, *] =f indgen (frame_no) 

;g=intarr (t_max) 

;q=[l, 0,0, 1,0, 0,1, 0,0, 1,0, 0,1] 

; correlation_state=l 

yes_no_values2=intarr (cell_no, f rame_no) 

for i=0, frame_no -1 do begin 

for j=0, cell_no -1 do begin 

yes_no__values2 [j , i] =-1 

endfor 

endfor 

for i=0, frame_no -1 do begin 
OK=0 

Sl for j=0, w_s -1 do begin ;possibly change the j range to make the window 
smaller 

W if(((i-j) ge 0) and ((i-j) It t_max) ) then begin 
0 if ( q(i-j) eq correlation__state ) then 0K=1 
Q endif 
fl| endfor 

\| for k=0, cell_no -1 do begin 

P if (OK eq 1) then Yes_no_values2 [k, i] =yes__no_values [k, i] 
*p. endfor 
endfor 

jj"* /print, 'frmae number is ' , frame_no 
f* ; print, 't max is ' , t_max 
tJ /print, 1 w_s is f , w_s 

m 

|,|, xx=where (yes_no_values2 [0 , *] eq -1, count) 

yes_ no_values3=intarr (cell_jno, f rame_jno - count) 

for j=0, cell_no -1 do begin 

k=0 

for i=0, frame_no -1 do begin 

if (yes__no__values2 [ j , i] ne -1) then begin 

yes_no_values3 [ j , k] = yes_no_values2 [ j , i] 

k=k+l 

endif 

endfor 

endfor 

/print, 'these are yes no values' 
/print, yes_no_values 
/print, q 

/print, 'these are yes no values 2' 
/print, yes_no_values2 
/print, 'these are yes no values 3' 
/ print , yes_no__values3 

print, 'State constrained yes and no values created' 
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draw_cross2 
end 

endcase 
end 

pro par_reestimation 

common mother__com, pixel_array, yes_no_values , coef , location, cell_no, 
frame_no, file_name, direct, $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common markov, 

n t m, A,B, P f observation, f , g, d, Fi, Q, back,xsi , Gamma, n_i , a__i , b_i,p_i , c , w_s , t_max 
t max=frame no - w s +1 



!except=2 

Let ' s recalculate forward variables : F 

f =dblarr (t_max,n) 

for i = 0, n-l do begin 

f (0,i) =P (i) *B (i, observation (0) ) 

endfor 

for t=0, (t_max-2) do begin 

for j=0, (n-l) do begin ; f (t , i) , A(i , j ) 

f (t+1, j) =sum(0,n-l, t , j , f , A) *B (j , observation (t+1) ) ;i goes from 0 to n-l 

endfor 

endfor 



I will apply scaling procedure to prevent underflows and increase accuracy 
c's will be the scaling coefficients 
f c = db 1 ar r ( t__max } 

f_init-dblarr (t_max,n) 
i A: f=dblarr (t_max,n) 
p s=double(0) 
ffj for i=0, n-l do begin 

p f_init (0 , i) =P (i) *B_prob (i , observation, 0 ,B, cell_no) 
|4 ; s=s+f_init (0, i) 

; print, 'state is f , i, 'b_prob is , B_prob (i , observation, 0 , B, eel l_no) 

endfor 

c [0]=l/s 

for i=0, n-l do begin 
f (0,i)=f_init(O,i)*c[0] 
endfor 



for t=0, (t_max-2) do begin 
s=double (0) 

for j=0, (n-l) do begin ; f (t , i) , A (i , j ) 

f_init (t+1, j ) =sum(0 , n-l, t , j , f , A) *B_prob ( j , observation, (t+1) , B, cell_no) ; i goes 

from 0 to n-l 

s=s+f _init (t+l, j ) 

endfor 

c [t+1] =l/s 

for j=0, (n-l) do begin 

f (t+1, j)=f_init (t+1, j) *c [t+1] 

endfor 

endfor 



NY02:354955.l 



-109- 



/Let's define the backward variables: BACK 
Back_init=dblarr (t_max,n) 
Back=dblarr (t__max,n) 
for 2=0, n-1 do begin 
Back_init (t_max -l,i)=l 

Back(tjmax - 1 , i) =Back_init (t_max -1, i) *c [t_jnax-l] 
endfor 

for t=t__max-2, 0, -1 do begin 
for 1=0, (n-1) do begin 
s=double (0) 
for j=0,n-l do begin 

s=s + A(i,j) *B_prob(j , observation, (t+1) ,B,cell_no) *Back(t+l,j) 
endfor 

Back_init (t , i) =s 

Back (t, i) =Back_init (t , i) *c [t] 

endfor 

endfor 

y ; ; print, 'these are backward variables' 
/print, Back 
/print, 'space' 

111 /Now let's define the sequence of probabilities xsi(t,i,j) - it is the 

"M /probaility of being at state i at time t and at state j at time t+1 

ffj xsi=dblarr (t_max, n, n) 

fn for t=0,t_max-2 do begin 

^ s=double(0) 

| ; .^ for i=0,n-l do begin 

I, for j=0,n-l do begin 

P s=s+F(t, i) *A(i, j ) *B_prob (j observation, (t+1) ,B, cell_no) *Back (t+1, j ) 

endfor 
IU endfor 

t«| for i=0,n-l do begin 
bk for j=0,n-l do begin 

xsi [t,i, j] = (F(t, i) *A(i, j) *B_prob(j observation, (t+1) ,b,cell_no) *Back(t+l, j) ) /s 

endfor 

endfor 

endfor 

/print, 'This is xsi' 
/print, xsi 
/print, 'space' 



/Gamma (t,i) is the probability of being at state i at t given the 

/observation and the model 

Gamma=Dblarr ( t_jnax, n) 

for t=0 , t__max -1 do begin 

for i=0, n-1 do begin 

s=double (0) 

for j=0, n-1 do begin 
s=s+xsi (t, i, j) 
endfor 

Gamma [t , i] =s 
endfor 
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endfor 

/print, 'This is gamma' 

; print, gamma 

; print, 'space' 



;Now we will begin the reestimation procedures of initial parameters P,A,B 



;reestimation of P 
for i=0, n-1 do begin 
P ( i ) =gamma ( 1 , i ) 
endfor 



/Reestimation of A 
for i=0, n-1 do begin 
for j=0, n-1 do begin 
sl=double (0) 
s2=double (0) 

for t=0, t__max-2 do begin 

sl=sl+Xsi (t, i, j ) 

s2=s2+Gamma (t, i) 

endfor 
M ; A(i, j)=sl/s2 
Q endfor 
p endfor 
m 



V: \ 



rii 



Reestimation of B 
for j=0, n-1 do begin 
for k=0, m-1 do begin 
sl^double (0) 
s2=double (0) 

for t=0, t_max-l do begin 

if (observation (t) eq k) then begin 

sl=sl+Gamma (t , j ) 

endif 

s2 = s2+Gamma (t , j ) 
endfor 

B(j f k) =sl/s2 

endfor 

endfor 

/Reestimation of B 
for j=0, n-1 do begin 
for k=0, m-1 do begin 
sl=double (0) 
s2=double (0) 

for t=0, t_max-l do begin 

sl=sl + Gamma (t, j )* (observation [k, t] ) 



s2=s2+Gamma (t, j ) 
endfor 

B(j ,k)=sl/s2 

endfor 

endfor 
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end 



; to sum n elements, used by hidden_markov_event . pro 



function sum, x,y,t,j,f,A 

; s will be used to denote the sum 

s =f (t,0)*A(0, j) 
for i=(x+l), y do begin 
s=s + f (t / i)*A(i / j) 
endf or 



return, s 
end 
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APPENDIX C 



; Name : medsub . pro 

/Description: This uses the running median filter to subtract the baseline 
;etc . 

/Variables: medsize is the size of the window used in calculating the filter 
; cellsin is a cells X frames array that contains the unfiltered 

data 

; cellsout has the cells after the base has been subtracted from 

them. 



pro medsub, cellsin, medsize, cellsout 

cellsout=f loat (cellsin) 
ncells=n_elements {cellsin (*, 0) ) 
nf rames=n_elements (cellsin (0, *) ) 
tmpl=f ltarr (nf rames+2*medsize) 
for i=0 , ncells-1 do begin 

tmpl (medsize :medsize+nf rames-1) =reform (cellsin (i, *) ) 
; tmpl= (tmpl-median (tmpl , medsize) ) >0 
tmpl=abs (tmpl-median (tmpl, medsize) ) 
jU, cellsout (i, *) =tmpl (medsize :medsize+nf rames-1) 

PI endf or 
M end 

pro nonnegf ac_converge, x, wl,vl 

Iteration for computing nonnegative factorisation 
x = Matrix; wl, vl = nonnegative factors. 



x = N x M 

Wl = N X R 
vl = R x M 



wv=wl#vl 



m xwv=x*0 

nz=where (x gt 0) 
J ;;; f' eps=le-20 

xwv(nz) =x(nz) / (wv(nz)+eps) 

vl=vl* (transpose (wl) #xwv) 



xwv=xwv*0 
wv=wl#vl 

xwv (nz) =x (nz) / (wv (nz) +eps) 



wl=wl* (xwv#transpose (vl) ) 

norm=rebin (transpose (total (wl, 1) ) +eps , n_elements (wl (* , 0) ) ,n__elements (wl (0 , *) ) ) 
wl=wl/norm 



end 



pro nonnegsvd__converge, mat, frac, u, v, nkeep, niter, errors 

;svd, mat, w, u, v 

; index=reverse (sort (w) ) 

; w=w (index) 

;u=u(*, index) 
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; v=v (* , index) 
sz=size (mat) 
x=sz (1) 
y=sz (2) 
R=nkeep 

noise=randomu (seed,R,y) 
noise2=randomu (seed, x, R) 
w=f Itarr (x,R) +. 01*noise2 
for i=0, R-l do w(i,i)=l 
v=mat (0:R-1, *) +. l*noise 

;U=U (*, 0:R-1) 

;V=v(*,0:R-l) 

;for t=0, R-l do begin 

;u(*,t)=u(*,t) *sqrt (w(t) ) 
;v(*,t)=v(*,t)*sqrt (w(t) ) 

; endf or 

; w=abs (u) 

; v=abs (transpose (v) ) 
errors=f ltarr (niter+1) 



for j = 0, niter do begin 
\& nonnegf ac_converge , mat , w,v 
fpl recon=w#v 
fii error =f ltarr (R) 
f% norm= total (mat A 2) 

p f or i=0,R-l do error (i) = total ( (mat-w (*, i) #v (i , *) ) A 2 ) /norm 
; 5? t! srt=sort (error) 
^ f rac=l -error (srt) 
W u=w (* , srt) 
IB v=v(srt, *) 

$ errors ( j ) = total ( (mat-recon) A 2) /norm 
\a endfor 

end 

M 

Ijl ;NAME : scalculate_svd 

IZ ; DESCRIPTION : using this, the user finds the nonnegsvd by inputing the number of 
ff iterations and 

;the number of modes to keep. (Though it is called calculate_svd / it is a 
different linear 

/factorization- nonnegative one.., but close to the general SVD) . 
pro scalculate_svd 

common svd__com, pixel_array, yes_no_values , coef, location, cell_no, frame_no, 
f ile__name , direct , $ 

total_f rame_no, time_re solution, x_size, y_size, box_size 
common pp_com, tp_pixel__array 

base=WIDGET_BASE (/column, title== 1 Perform Non-negative SVD') 
draw=WIDGET_DRAW (base , xs ize=4 0 0, ys i ze=4 0 0 ) 

modeskeep=CW_FIELD (base, title= 1 Enter the number of modes that should be 
kept : 1 , $ 

VALUE^cell^no, /integer, /return_events) 
niterations=CW_FIELD (base, title= ' Number of iterations:', $ 

VALUE=100 , /integer, /return_events) 
button=WIDGET_BUTTON (base, value= 1 Calculate SVD ' , uvalue= ' button_pressed ' ) 
widget_control, /realize, base 

state={modeskeep :modeskeep, niterations : niterations } 
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widget_control, widget_info (base, /child) , set_uvalue=state 

xmanager, ' scalculate_svd' , base 

end 

;NAME : scalculate_svd__event 

/DESCRIPTION: This is the event handler for scalculate_svd . . It calls 
nonnegsvd_converge which calls nonnegf ac_converge based on the 

;number of iterations and the number of modes to keep. **Make sure to do all the 
required error checking...** 

pro scalculate_svd_event , event 

common pp_com, tp_pixel_array 

common svds__com, w, u, v, errors, mkeep, imu 

common svd_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, 
f ile_name, direct, $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
common flags, svd_calculated, cell s__de fined 
common vis_svdcom, visbase, visbase2 

common local, imvar ; unused 

commons such as this should be removed ;this makes it readable from 'Choose 
Cells' too 

common stockdata, symbol_array, date_array 



m 



if this is 1, the general 
if this is 1, the spatial 
find what is in /Child 



stateholder=WIDGET_INFO (event .handler , /Child) 
print, visbase, visbase2 

vismanaged=VJIDGET_INFO (visbase, / MANAGED ) 
visualization widget is open 

vismanaged2=WIDGET_INF0 (visbase2 , /MANAGED) 
visualization widget is open 

WIDGET_CONTROL, stateholder, GET_UVALUE= state 
WiDGET_CONTROL , state . modeskeep , GET_VALUE=mkeep 
WIDGET_CONTROL , state . niterat ions , GET_VALUE=niter 
Widget_control, event . id, get_uvalue=uval 

;if ( (vismanaged eq 1) and (mkeep ne mkeepold) ) OR ( (vismanaged2 eq 1) and 
(mkeep ne mkeepold2)) then begin ; number of modes being calculated is 

not number being showed in visualization. general 
4 ; mess=WIDGET_MESSAGE ( ' A visualization widget is open with a capacity 

'J for a different number of modes than you currently specify] Change the number 
of modes, or close the visualization widget and repeat.', /ERROR) 
if ( (vismanaged eq 1) OR (vismanaged2 eq 1) ) then begin 

mess=WIDGET_MESSAGE ( 1 To avoid conflicts with the number of modes 
specified, please shut your visualization widgets . ' , /INFORMATION) 
endif else begin 

if uval eq ' button_j?ressed ' then begin 
if (mkeep gt cell_no) then begin 

me ss=D I ALOGJVIES SAGE ( 'Number of modes cannot be more than the number 

Of Cells! ' , /ERROR) 

endif else begin 

nonnegsvd_converge, tp_pixel_array, w, u, v, mkeep, niter, errors 
; imu=f ltarr (x_size,y_size, mkeep) 
;for i=0, mkeep- 1 do begin 

; f or this to work with stocks -> 

; iml=f ltarr (x_size , y_size) 

; imu (* , * , i) =imvar 

; iml (location (*) . coord (0) , location (*) . coord (l) ) =u (*, i) 
; reconstruct to dimensions of image 

;for j=0, cell_no-l do begin ;draw the 

boxes 
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; imu (location (j ) . coord (0)- 

box_size: location (j ) . coord (0) +box_size, location (j ) . coord (1) - 

box__size : location (j ) .coord (1) +box_size, i) =iml (location (j ) . coord (0) , location (j ) . c 
oord(l) ) 

; endf or 

/endf or 

plot, w, xtitle='Mode Number', yt it le= ' Power Contribution' 
svd_calculated-l 
endelse 
endif 

endelse 
end 



; NAME: sload_and_convert_excelf ile 

; ; This file reads MS Excel file with stock data of the following format 
; ; date tickerl ticker2 ... ... 

; ; value closel close2 ... ... 

;; This file creates the following arrays: 

; ; symbol_array (string) - first row of the excle file minus first value 
Lj,; ; date_array (long or string or date) - first column of the excel file minus 
f|the first value 

p; ; pixel_array (float) - all the rest 

J£; ; pixel_array later gets transformed by calculating deltaF/F - under the same 

rS name 

Nipro sload_and_convert_excelf ile 

COcommon svd_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, 
■Lfile_name, direct, $ 

% total_jErame_no, time_resolution, x_size, y_size, box_size 

common mother_com, pixel_array, yes_no_values , coef, location, cell_no, 
Irframe no, file name,$ 

total_f rame_no, time_resolution, x_size, y_size, box_size 
^/common test, s trying, state3 
|jcommon old_skool_data, original__data 

common flags, cells_def ined, spikes_def ined, single_jplot_def ined, 
Hcor r e l_coe f _de f ined 
common flags, svd_calculated, cells_def ined 

common with_choose_cells__com, thresholdl, threshold2 , rms_threshold 
common textf ile_vars , text_flag, diode_array, max_jium_of_diodes 
common stockdata, symbol_array, date_array 



mes sage_dialog=WIDGET_MES SAGE ( "This procedure loads stock data, it creates 
symbol_array, date_array, and pixel_array" , / INFORMATION) 

f ilename=dialog_j?ickf ile (/read, f ile= ( 1 stocks . slk' ) , get_path=f ilepathl) 
if (filename eq 1 1 ) then begin 

message_dialog==WIDGET_MESSAGE ("No data read.", /INFORMATION) 
endif else begin 

close, /all 

openr, 1, filename 



/use excel import function to fill all the required arrays 
symbol_array= read_sylk (filename, / ARRAY , nrows = 1 , startcol=l) 
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date_array=read_sylk (filename, /ARRAY, 110018 = 1, 
start row=l, /uselongs) 

data = read_sylk (filename, / ARRAY, startrow=l, startcol=l) 
close, /all 

cell_no = n_elements (symbol_array) 
frame_no = n_elements (date_array) 
; state3={ frame : frame} 

;WIDGET_CONTROL, state3 .name, SET_VALUE=strmid (filename , 
strlen(f ilepathl) ) 

; WIDGET_CONTROL , s tate3 . frame , SET_VALUE=f rame_no 
;WIDGET_CONTROL, state3 .time, SET_VALUE=time_re solution 

;; build 'pixel_array' 

pixel_array = f ltarr (cell_no, frame_no) 
for j=0, frame_no -1 do begin 

for i = 0, cell_no - 1 do begin 

pixel_array [i, j] = data[j,i] 

endfor 

endfor 

original_data = pixel_array 
;the following finds delta F over F 
deltaf =f ltarr (cell_no, f ramejio) 

for i = 0, cell_no - 1 do begin 

for j=l, frame_no-l do begin 

; deltaf [i, j] =100* (pixel_array [i , j ] - pixel_array [ 

array [i, j-1] 

deltaf [i, j] =pixel_array [i , j ] 

endfor 

endfor 

for i = 0, cell_no - 1 do begin 
deltaf [i,0] = 0 

endfor 

pixel_array=deltaf 

; / initialize other variables 
yes_no_values = intarr (cell_no, frame__no) 
coef = f ltarr (cell_no, cell_no) 
rms_threshold = f ltarr (cell_no) 
thresholdl = f ltarr (cell_no) 
threshold2 = f ltarr (cell__no) 
for cell = 0, cell_no - 1 do begin 
rms_threshold [cell] =2.0 
thresholdl [cell] =2.0 
threshold2 [cell] =3.0 

endfor 

cells_def ined = 1 
; spikes_def ined = 0 
; single _j?lot_de fined = 0 
; correl_coef_def ined = 0 
;box_size = 2 

;x_size = max (x_pos_data) + 350 
;y_size = max (y_jpos_data) + 3 50 
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;; locations are shifted up and to the right by 300 so that all 
locations are positive (and thus can be plotted correctly) 

/location = replicate ( {struct , coord : intarr (2) , size:0, 
half_side:0 . 00} , celljio) 

rfor i = 0, (cell_no - 1) do begin 
location (i) . size = cell_no 
location (i) . half_side = 1 

location (i) . coord [0] = xjpos_data [diode_array [i] - 1] + 300 
location (i) .coord [1] = y_j)os_data [diode_array [i] - 1] + 300 
; endf or 

text_flag = 1 

message_dialog=WIDGETJVIESSAGE( "Finished reading data.", 

/ INFORMATION) 
endelse 

end 

;Name: spreprocess_wid 

;Decription: This is the widget that allows you to preprocess your data 
(smooth, medsub, abs) 
/Modification History 

H ; 

£| pro spreproces s_wid 

p common svd_com, pixel_array, yes_no_yalues, coef, location, cell_no, 

ft% f r ame_no , f i 1 e_name , di r e c t , $ 

|p total_f rame_no, time_resolution, x_size, y__size, box_size 

!j f t base3=WIDGET_BASE {/column, title= 1 Preprocess ' ) 

,*f draw5=WIDGET_DRAW(base3,xsize=550, ysize=300) 

slide2=WIDGET_SLIDER (base3 , minimum=l / maximum=cell_no, title= ' Stock 
63 Number 1 ) 

ifi button=cw_bgroup (base3 , /row, [ 1 Plot 1 , 'Median Subtract', 'Smooth' , 'Take 

Absolute Value 1 , 'Return Original Values 1 ], /return_name, UVALUE=0) 
|4 ; stockwin=CW_FIELD (base3 , title= 1 Stock 1 ) 

m medwin=CW_FIELD (base3 , title= ' Running Window Width for Median 

i|f Subtraction ' , VALUE=2 0 , / INTEGER , / RETURN_EVENT S ) 

smoothwin=CW_FIELD (base3 , title=' Width of Smoothing Window 
,VALUE=3, /INTEGER, / RETURNEE VENTS ) 

statel= { slide : slide2 , medwin rmedwin, smoothwin : smoothwin, draw : drawS } 
widget_control, base3 , /realize 
widget_control , base3 , set__uvalue=statel 
xmanager, ' spreprocess_wid ' , base3 

end 

; Name : spreproces s_wid_event 
/Description: Widget handler 

pro spreproces s_wid_event , event 

common svd__com, pixel_array, yes_no__values, coef, location, cell_no, frame_no, 
f ile_name, direct, $ 

total_f rame_no, time__resolution, x_size, y_size, box_size 
common pp_com, tp_j?ixel_array 

widget_control , event. top, get_uvalue=state3 

WIDGET_C0NTR0L, state3 . slide, GET_VALUE=:cell_numb 

WIDGET_CONTROL, state3 . smoothwin, GET_VALUE= swin 

WIDGET_CONTROL, state3 .medwin, GET_VALUE=mwin 

WIDGET_C0NTR0L, state3.draw, GET VALUE = win num 



NY02:354 955.1 



-118- 



;WIDGET_CONTROL, state3 . stockwin, SET_VALUE=symbol_array (cell_numb-l) 
if (mwin le 1) then begin 

mess=WIDGET_MESSAGE ( 'Median filter window width must be greater than 

II', /ERROR) 

endif else begin 

if ( (swin It 2) OR (swin ge f rame_no) ) then begin 
print, frame_no 

mess =WIDGET_MESSAGE (' Smoothing window width must be >2 and <Total 
number of frames I ' , /ERROR) 
endif else begin 

if { (size (event .value) ) [1] eg 7) then begin ;; i.e. if we've pressed a 
button which returns a string as its value... 
Case event. value Of 



1 Smooth > : BEGIN 
wset, win_num 

for i=0, cell_no-l do begin 

tp _jpixel_array (i , *) = smooth (tp_pixel_array (i , *) , swin) 

endf or 

; p 1 o t , f indgen ( f r ame_no ) * t ime_r e so lu t i on , tp_p ixe l_ar r ay (eel l_numb - 

l f *) 

| r | ; plot, f indgen (frame_no) , tp_pixel_array (cell_numb-l f *) , 

C: 



l^ytitle=' Delta F/F 1 , xtitle= 1 Days ' 



1 :^ 



END 



■Median Subtract 1 : BEGIN 



¥1 ii 

wset, win__num 

h "4 medsub, tp_pixel_array , mwin, tp_pixel_array 

to 
m 

y plot, f indgen (frame_no) , tp_pixel_array (cell_numb-l , * ) , 

| 4 ytitle= 'Delta F/F', xtitle= ' Days 1 



;plot , f indgen (f rame_no) *time_resolution, tpj?ixel_array (cell_numb- 



p 
m 

H 



W wset, win num 



END 

'Take Absolute Value 1 : BEGIN 



tp_jpixel_array=abs (tp_pixel_array) 

;plot , f indgen (f rame_no) *time_resolution, tp_pixel__array (cell_numb- 



1,*) 

plot, f indgen (frame_no) , tp_pixel_array (cell_numb-l, *) , 
ytitle= 'Delta F/F', xtitle= ' Days ' 
END 



'Return Original Values': BEGIN 
wset, win^num 
tp_pixel_array=pixel_array 

/plot, f indgen (f rame_no) *time_resolution, tp_pixel_array (cell_numb- 

!,*} 

plot, f indgen (frame_no) , tp_pixel__array (cell_numb-l, *) , 
ytitle=' Delta F/F' , xtitle= 1 Days ' 
END 



■ Plot » : BEGIN 
wset , winjium 

;plot , f indgen (f rame__no) *time_resolution, tpjpixel_array (cell_numb- 

1,*) 
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plot, f indgen (frame_no) , tp_pixel_array (cell_numb~l, *) , 
ytitle=' Delta F/F 1 , xtitle= ' Days ' 
END 

ELSE : 
END CASE 
endif else begin 

wset, win_num 

;plot, f indgen (frame_no) *time_re solution, 
tp_j)ixel_array (cell^numb-l, *) 

plot, f indgen (f rame_no) , tpjpixel_array (cell_numb-l, *) , 
ytitle= 'Delta F/F', xtitle= 'Days 1 
endelse 
endelse 
endelse 
end 

; Name: s superimposed _jplots 

; Description: This program plots the deltaF/F values of a single cell per time. 

pro ssuperimposed__plots, delta__f_values , celljnumber, frame_number 

jp; !p .multi=0 
f! window, 2 

Q ;plot, f indgen (frame_number) *time_resolution, delta_f_values {0, *) , 

fiyrange= [min (delta_f_values) -2 , max (delta_f_values) +2] , $ 

j|t ; title='All cells', ytitle= 1 iEF/F ' , xtitle='time in seconds' 

l* 4 ! plot, f indgen (f rame__number) , delta_f_values (0 , *) , yrange= [min (delta_f_values) - 

2,max(delta_f_yalues) +2] , $ 
W title='All stocks', ytitle= ' JEF/F ' , xtitle=' Daily market closes' 

ISxyouts, -6, 0, 1 
3 if cell_number gt 1 then begin 
|4 for cell_numb=l , cell_number- 1 do begin 
|4 ; oplot , f indgen (frame_number) *time_resolution, 

(delta_f_values (cell_numb, *) ) 
■jijlj oplot, f indgen (f rame_number) , (delta_f_values (cell_numb, *) ) 

]% ;xyouts, -6, cell__numb*20, cell_numb+l 

a= (delta_f_values (cell_numb, *) eq max( (delta_f_yalues (cell_numb, *) ) ) ) 
^ a=where (a eq 1) 

xyouts, a,max( (delta__f_values (cell_numb, *) ) ) , cell_numb+l 

endf or 

endif 

end 

; NAME : SSVD_gui 

; DESCRIPTION : This widget will be used to load data, pre-process, 
/perform an SVD, and then visualize the modes created. 

pro ssvd_gui 

common svd__com, pixel_array, yes_no__values , coef, location, cell_no, frame_no, 
file_name, direct, $ 

total_f rame_no, time_re solution, x__size, y_size, box_size 
common flags, svd_calculated, eel ls_de fined 
common vis_svdcom, visbase, visbase2 

visbase=long (-1) 
visbase2=long ( -1) 

svd__calculated=0 ;set flags to 0 
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cells_def ined=0 

base=WIDGET_BASE (/ column, title= * SVD STOCK GUI') 
stuf f = { cw_jpdmenu_s , flags : 0 , name : 1 1 } 
details=[{ cw_pdmenu_s, 0, 'Previous Experiment 1 }, $ 
;{ cw_pdmenu_s , 0, 'Choose Cells' }, $ 

{ cw_jpdmenu_s , 0, 'Load Data'}, $ ;for the new import of stock data 
{ cw_j?dmenu_s , 0 , ' Superimpose Traces ' } , $ 
{ cw__pdmenu_s , 0 , ' Preprocess 1 } , $ 
{ cwjpdmenu_s , 0 , ' SVD 1 } , $ 
{ cwjpdmenu_s , 1, 1 Visualize '} , $ 
{ cw_jpdmenu_s / 0, 'Across Stocks'}, $ 
{ cw_pdmenu__s , 2, 'Within Modes '} ] 
; { cwj)dmenu_s, 0, 'General'}^ 
;{ cw_j?dmenu_s , 2, 'With Locations'}] 
pull_down=cw_pdmenu (base , details , / return_f ull_name , UVALUE=12 ) 

stuf f 2= {cw_pdmenu_s2 , flags : 0 , name : ' ' } 
details=[{ cw_j)dmenu_s , 0, 'Save Analysis'}, $ 
{ cw_pdmenu_s, 0, 'Color'}, $ 
{ cw_j>dmenu_s , 0, 'Retail' }, $ 
{ cw_jpdmenu_s , 0, 'Exit'}] 
j-si pull_down2=cw_jpdmenu (base, details, /return_full_name, UVALUE=13) 

M ! 

% ;frameid=CW_FIELD(base, title='Total Number of Frames ■ , VALUE=total_f rame_no, 

W / INTEGER , $ 

l J* ; /RETURN_E VENTS ) 

" ( i ; timeid=CW_FIELD (base, title='Time Resolution in seconds ' ,VALUE=time_resolut ion, 

ffl /FLOATING, / RETURN_E VENTS ) 

C0 /don't need time resolution for this.. 

$ / state={ frame : frameid, timeitimeid} 

widget_control, /realize, base 
y ; ;widget_control, widget_inf o (base, /Child), SETJJVALUE=state 

widget_control , widget_inf o (base, /Child) 
p! xmanager, 1 ssvd_gui ' , base 
^ end 

fv% - — _ . 

; Name : s s vd_gui_e vent 

^* ; Synopsis: ssvd__gui_event , event 

; Description: This is the event handler for the ssvd_gui window. 

pro ssvd_gui__event , event 

common svd_com, pixel_array , yes_no_values , coef, location, cell_no, frame_no, 
file_name, direct, $ 

total_f rame_no, time_resolution, x_size, y__size, box_size 

common pp_com, tp_pixel_array 

common flags, svd__calculated, cells__def ined 

common stockdata, symbol_array, date_array 

stateholder=WIDGET_INFO (event . handler , /Child) 

; WIDGET_CONTROL , stateholder, GET_UVALUE = s t a t e ;find what is in /Child 
; WIDGET_CONTROL , state . frame , GET_VALUE=total_f rame_no 
;WIDGET_CONTROL, st ate. time, GET_VALUE=time_resolution 
Case event. value OF 

'Previous Experiment ': BEGIN 
; tempstate=state 
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data_file=dialog__pickf ile (/read, title= 1 Select the saved variables 
file', GET_PATH=filepathl) 

if (data_f ile eq 1 1 ) then begin 

mess=WIDGET_MESSAGE ( 'No saved variables file specified 1 , 

/ INFORMATION) 

endif else begin 

restore, data_file 

mess=WIDGET__MESSAGE ( 'Done loading variables', / INFORMATION) 
cells_def ined=l 

;WIDGET_CONTROL, temps tat e . frame, S E T_VALUE = t o t a 1 _f r ame_no 
; WIDGET_CONTROL , tempstate . time, SET_VALUE=time_resolution 
;not needed anymore for stocks 
endelse 

END 

'Load Data' : BEGIN 

sload_and_convert_excelf ile 

END 

' Superimpose Traces 1 : BEGIN 

;if { (total_frame_no eq 0) or (time_re solution eq 0)) then begin 

; mess=WIDGET_MESSAGE ( ' You must enter the total number of frames 

fli and time resolution before proceeding I f , /ERROR) 

|'4 ; endif else begin 

^ if (cells_def ined eq 1) then begin 

5 * tp _p ixe l_ar r ay =p ixe l_ar r ay 

; superimposedjlots, tpjpixel_array , cell_no, frame_jno, 

time_re solution 

ttj ssuperimposed_plots , tp__pixel_array, cell_no, frame_no 

|0 endif else begin 

$ me ss=WIDGET_MESS AGE ( 'Buddy, load data before plotting.', 

jri /ERROR) 

p 

J lfi f ' Preprocess 1 : BEGIN 

; if ( (total_f rame__no eq 0) or (time_resolution eq 0)) then begin 

mess=WIDGET__MESSAGE ( 1 You must enter the total number of frames 
and time resolution before proceeding!', /ERROR) 
; endif else begin 

if <cells_def ined eq 1) then begin 

spreprocess_wid 
endif else begin 

mess=WIDGET_MESSAGE ( 'Load data before plotting.', /ERROR) 
; endelse 
endelse 

END 

' SVD ' : BEGIN 

/ if ( (total_f rame_no eq 0) or (time_resolution eq 0)) then begin 
/ mess=WIDGET_MESSAGE ( ' You must enter the total number of frames 

and time resolution before proceeding!', /ERROR) 
endif else begin 

if (cells_def ined eq 1) then begin 

scalculate_svd 
endif else begin 



; endelse 
endelse 
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mess=WIDGETJVIESSAGE ( ' Load data before calculating the SVD.', 

/ERROR) 

;endelse 
endelse 

END 

; 'Visualize .General 1 : BEGIN 
r Visualize .Across Stocks': BEGIN 

;if ( (total_f rame_no eq 0) or (time_resolution eq 0)) then begin 
/ mess=WIDGET__MESSAGE ( ' You must enter the total number of frames 

and time resolution before proceeding! 1 , /ERROR) 
;endif else begin 

if (svd_calculated eq 1) then begin 

svisualize_wid 
endif else begin 

mess=WIDGET_MESSAGE ( 1 You must calculate the SVD before 
visualizing the modes . 1 , /ERROR) 
; endelse 
endelse 

END 

'Visualize. Within Modes': BEGIN 

if (svd_calculated eq 1) then begin 

svisualize_wid_locs 
endif else begin 

me s s =WIDGET MESSAGE ( ' You must calculate the SVD before 



visualizing the modes. ' , /ERROR) 
'' f 4 endelse 
CO END 

Co 

u 
u 

0 

^ data_file=pickfile (/write, f ile= ' Variables_svd#_' , title= 1 Create the 

H : saved variables file') 

if (data_file eq ' 1 ) then begin 

mess=WIDGET_MESSAGE ( 'No saved variables file specified 1 , 



'Color' : xloadct 
•Retail 1 : retail 
'Save Analysis': BEGIN 



/INFORMATION) 

/INFORMATION) 
END 



endif else begin 

SAVE, /VARIABLES, FILENAME=data_f ile , all, /verbose 
mess=WIDGET_MESSAGE ( 'Data from this experiment saved' , 

endelse 



'Exit 1 : BEGIN 

WIDGET_CONTROL, /DESTROY, event . top 

END 



END CASE 

end 



; Name: svisualize wid 
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; Description: This is the widget that shows the spatial and temporal modes of 
the data. 

pro svisualize_wid 

common svd_com, pixel_array, yes_no_valu.es , coef, location, cell_no, frame_no, 
file_name, direct, $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 

common svds_com, w, u, v, errors, mkeep, imu 

common wid, num_jplots 

common vis_svdcom, visbase, visbase2 

num_plots=0 ;no plots have been done, so the first must 

be a plot, not Oplot 

visbase=WIDGET_BASE (/column, title= 'Visualization Widget r ) 

drawl =WIDGET_D RAW (visbase , xsize=500 , ysize=10 0 ) 

draw2=WIDGET_DRAW (visbase, xsize=500 ,ysize=150) 

slide=WIDGET_SLIDER (visbase, value=l, maximum=mkeep , minimum=l, $ 

title- 1 Mode Number ' , uvalue= 1 slider__pressed ' ) 
buttonl=WIDGET_BUTTON (visbase, value='Plot' , uvalue= 'plot ^pressed 1 ) 
draw3 =WIDGET__DRAW ( vi sbase , xsi ze=5 0 0 , ysi ze=15 0 ) 

slide2=WIDGET_SLIDER (visbase, value=l, maximum=cell_jio, minimum=l,$ 

title=' Stock Number' , uvalue= ' slider2_pressed' ) 
butt on2=WIDGET_BUTT0N(visbase,value=' Add Stock to Plot uvalue =' Add 1 ) 
button3=WIDGET_BUTT0N (visbase , value= ' Clear Stocks ' , uvalue= f Clear ' ) 
WIDGET_CONTROL, /realize, visbase 
loadct, 2 7 

statel= {drawl : drawl , draw2 : draw2 , draw3 : draw3 , slide : slide, slide2 : slide2 , 
but ton2 : butt on2 , $ 

button3 :button3 } 
WIDGETjCONTROL , visbase , set_uvalue=statel 
xmanager, 1 svisualize_wid ' , visbase 

end 

; Name: svisualize_wid_event 

/ Description: This reads the mode chosen and displays the spatial and temporal 
modes for it. 

pro svisualize_wid_event , event 

common svds_com, w, u, v, errors, mkeep, imu 

common wid, num_plots 

common svd_com, pixel_array, yes_no_values , coef, location, cell_no, frame^no, 
file_name, direct, $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 

common pp_com, tp_j>ixel_array 

common stockdata, symbol_array , date_array 

WIDGET__CONTROL , event . top , get_uvalue= state 
WIDGET_CONTROL , state . drawl , get_value=winl 
WIDGET_CONTROL , state . draw2 , get_yalue=win2 
WIDGET_CONTROL , state . draw3 , get__value=win3 
WIDGET_CONTROL , state . slide, get_value=mode_no 
WIDGET_CONTROL , state . slide2 , get_value=p lot cell 
WIDGET_CONTROL, event . id, get_uvalue=uval 

]y.minor=-l ;this is so that you can see the values on the y axis- 
suppress the tick-marks 

if ( (string (uval) eq 'plot ^pressed' ) OR (string (uval) eq 
1 slider ^pressed' ) ) then begin 

wset, win2 & plot, v (mode_no-l , *) , xtitle= ' Days ' , 
ytitle=' Amplitude 1 , t it le= r Temporal Mode' 
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wset, winl & plot, f indgen (cell_no) +1 , u (* , mode_no-l) , psym=2, 
xtitle=' Stock Number', title= ' Spatial Mode' 

;wset, winl & plot, u(*,mode_no-l) , 
Xtickname=syrnbol_array,xticks=cell_no, psym=2, xtitle= ' Stock ' , title= ' Spatial 
Mode ' 

;wset, winl & bar__plot, u (* , mode__no-l) , barnames=symbol_array 
num_plots=0 

endif 

if (string (uval) eg 'Add') then begin ;Add 
Cell to Plot was pressed 

if (num_plots eg 0) then begin 

wset, win3 & plot, tp__pixel_array (plotcell-l , *) , 
title= l Superimposed Plots' 

num_plots=l 
endif else begin 

wset, win3 & oplot, tp_pixel_array {plot cell -1, *} , 
color=num_plo t s * 1 0 0 

n um_jp 1 o t s = num_j) 1 o t s + 1 
endelse 

endif 

y ; if (string (uval) eg 'Clear') then begin /Clear Cells 

m was pressed 

pi num__plots = 0 

wset, win3 & erase 

^ endif 

! ^ iy.minor=0 

'"4 end 

co 

l9 ; Name: svisualize_wid__locs 

3i ; Description: This shows the spatial modes of the components with specific 
|4, locations 

f% pro svisualize_wid_locs 

pll common svd_com, pixel_array, yes_no_values , coef, location, cell_no, frame_no, 
file_name, direct, $ 

total_f rame_no, time_re solution, x_size, y_size, box_jsize 
^ common svds_com, w, u, v, errors, mkeep, imu 
common vis_svdcom, visbase, visbase2 
common Iocs, n_rows, n_columns 

n_columns=round(sqrt (cell_no) ) ; numb e r of columns to make 

closest to sqre 

n_rows=ceil (cell_no*l. 0/n__columns) 

visbase2=WIDGET__BASE (/column, title= ' Scroll through the spatial modes') 
;draw=WIDGETJDRAW(visbase2, xsize=x_size , ysize=y_size) 
draw=WIDGET_DRAW(visbase2, xsize=n_rows*7 0, ysize=n_columns*7 0) 
slide=WIDGET_SLIDER(visbase2, value=l, maximum=mkeep, minimum=l, 
title='Mode Number', $ 

uvalue= 1 slider_j>ressed f ) 
buttonl=WIDGET_BUTTON(visbase2, values' Plot 1 , uvalue= 1 plot_pressed 1 ) 
WIDGET__CONTROL , /realize, visbase2 
statel={ draw: draw, slide : slide} 
WIDGET_CONTROL, visbase2 , set__uvalue=statel 
loadct, 5 

xmanager , ' svisualize_wid_locs ' , visbase2 

end 
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; Name: svisualize_wid_locs_event 

/ Description: This will show the spatial locations of the cells from the modes 
pro svisualize__wid_locs_event , event 

common svd_com, pixel_array, yes_no_values , coef, location, cell_no, frame_no, 
f ile_name, direct, $ 

total_f rame_no, time_resolution, x_size, y_size, box_size 

common stockdata, symbol_array , date__array 

common svds_com, w, u, v, errors, mkeep, imu 

common Iocs, n_rows, n_columns 

cell_numb=0 

WIDGET_COMTROL, event . top, get_uvalue=state 
WIDGET_CONTROL , state . draw, get_value=win 
WIDGET_CONTROL., state, slide, get_value=mode__no 
WIDGET^CONTROLi, event, id, get_uvalue=uval 
if ( (string (uval) eq ! plot_pressed 1 ) OR (string (uval) eq 
' slider_pressed' ) ) then begin 
;loadct, 5 
wset, win 

;the best way to represent the data seems to be in a matrix built of 

symbols 

^ /with round (sqrt (totalstocks) ) as the number of columns 

CI /tvscl, imu(*, *,mode_no-l) 

Cl imu=f ltarr (2 00*n_columns , 250*n_rows) 

indexl= (10+n_columns-l) *n_columns 

;for i=20, 5*indexl, 5*n__columns do begin ; to give it space 

uto draw a box 

f$ ; for j=20, (10+n_rows-l) *5*n_rows, 5*nj:ows do begin 

l^j ;for i=10, (10+n_columns-l) *10, 10 do begin 

* w for i=50, (l+n_rows-l) *50, 50 do begin 

f for j=60, (l+n_columns-l) *60, 60 do begin 

* eel l_numb = c e 1 l_numb + 1 

I"* if (cell_numb le cell_no) then begin 

Q imu (i-5 : i+5 , j -5 : j +5) =u (cell_numb- 1 , mode^no- 1 ) ; specify 

Jljbox size later 

|3i end if 

|4, ;xyouts, i-50, j-25, cell_numb, /device 

endfor 

endfor 

/for i=0, celljno-1 do begin 

xyouts, location (i) . coord (0) , location (i) . coord (1) , 
string (i+1) , color=100 , /device 
; endfor 

tvscl, imu /now write stock symbols on top 

cell_numb=0 

for i=50, (10+n_rows-l) *50, 50 do begin 

for 3-60, (l+n_columns-l) *60, 60 do begin 

ce 1 l_numb=ce 1 l_numb+ 1 

/print, cell_numb 

if (cell_numb le cell_no) then begin 
/xyouts, i-50, j-25, cell_numb, /device 
xyouts, i-35, j+15, cell_numb, /device, 

charsize= . 8 

xyouts, i-15, j-25, symbol_array (cell_numb-l) , 
/device, charsize=1.2 / , color=u (cell_numb-l , mode_no-l) *5 
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endif 
endfor 

endf or 

endif 



end 



ru 

IB 
CO 
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